We study the life-cycle of miscible fingering, from the early fingering initiation, through their growth and nonlinear interactions to their decay to a single finger at late times. Dimensionless analysis is used to relate the number of fingers, the nature of their nonlinear interactions (spreading, coalescence, tip splitting), and their eventual decay to the viscosity ratio, transverse Peclet number, and anisotropic dispersion. We show that the initial number of fingers that grow is approximately half that predicted by analytical solutions that neglect the impact of longitudinal diffusion smearing the interface between the injected solvent and the displaced fluid. The growth rates of these fingers are also approximately one quarter that predicted by these analyses. Nonetheless, we find that the dynamics of finger interactions over time can be scaled using the most dangerous wavenumber and associated growth rate determined from linear stability analysis. This subsequently allows us to provide a relationship that can be used to estimate when predict when the late time, single finger regime will occur.

Viscous fingering, also known as the Saffman-Taylor instability (Ref. 1), occurs when a less viscous fluid displaces a more viscous fluid in a porous medium. It can be observed in many applications including chromatography (Ref. 2), geological carbon dioxide (CO2) sequestration (Ref. 3), and enhanced oil recovery processes (Ref. 4). In chromatography and oil recovery, the aim is to reduce or prevent its occurrence as it reduces separation efficiency (chromatography) and oil recovery; however, fingering may be beneficial in subsurface geological storage of carbon dioxide by increasing mixing and hence trapping of the carbon dioxide (Refs. 5 and 6).

It is challenging to model viscous fingering numerically, especially in miscible displacements, as the pattern and growth of the fingers depend on physical diffusion and dispersion as well as the viscosity ratio between the displacing and displaced fluids (Refs. 7–10). The larger the transverse dispersion, the smaller the number of fingers that grow initially (Refs. 8 and 11). Simulations typically require very fine meshes and/or higher order methods to ensure physical diffusion and dispersion dominate over truncation errors and predict the correct number of fingers and their growth rates (Refs. 12–14). As a result, many authors have investigated viscous fingering as a means for demonstrating the utility of their numerical schemes (Refs. 15–17).

The dynamics of viscous fingers are highly nonlinear and cannot be described analytically except at early time when the fingers grow independently of their neighbors. In this case, the growth rate and number of fingers can be estimated using perturbation analysis (e.g., Refs. 8, 11, 18–20, and 21). At intermediate times, numerical simulations and/or physical experiments are required to determine the complex interactions between the fingers that occur as they merge and split in a nonlinear fashion until at late times only a single finger remains (Refs. 11, 21, and 22).

Although there is a wide literature on viscous fingering in porous media, most investigations have focused on the early time regime (when the system is amenable to mathematical analysis). Some studies have investigated fingering at intermediate times (when there are multiple fingers interacting), but few have quantified the change in the number of fingers and identified the mechanisms that operate during the decay, merging, and splitting stages as a function of the viscosity ratio and dispersion. Only Tan and Homsy,7 Zimmerman and Homsy,11,21 and Nijjer et al.22 have studied the full life cycle of viscous fingering, and these studies considered systems with periodic boundary conditions perpendicular to the average flow direction. Morris and Ball23 and subsequently Yang and Yortsos24 have both shown that fingering patterns may be altered by the presence of no-normal flow boundaries, especially for aspect ratios greater than 1. In particular, Yang and Yortsos24 observed that any fingers growing adjacent to the boundary may grow more quickly than fingers in the interior of the system. Most real systems (chromatography columns, aquifers being used for carbon dioxide sequestration and oil reservoirs) have no-normal flow boundaries.

The average behavior of the flow in the different fingering regimes and the times at which the different regimes start are of particular interest to engineers as they typically use empirical models such as those of Koval25 or Todd and Longstaff26 to predict the impact of fingering on the trapping of carbon dioxide or oil recovery. These empirical models assume that the average concentration profile of the displacing fluid grows linearly with time. Various studies (Refs. 27–29) have shown that these empirical models provide a good match to both experimental data and detailed simulations of viscous fingering at modest viscosity ratios. There is, however, some evidence that the average concentration of the displacing fluid is not self-similar with respect to distance/time, as assumed by these models, for larger mobility ratios (Ref. 30), for low Peclet numbers (Ref. 31), and in the late time, single finger regime (Ref. 32).

In this paper, we study the different fingering regimes observed in unstable, miscible displacement in 2D, homogeneous, rectilinear systems with no-flow boundaries perpendicular to the average direction of flow. The impacts of mobility ratios M, transverse diffusion DT, and longitudinal diffusion DL on the flow regimes are investigated. It is important to investigate simplified 2D systems before investigating more realistic 3D systems in order to ensure that we have a full description of the appropriate physics in these systems. 3D systems are both more complicated to analyze mathematically and more numerically intensive to simulate; moreover, previous studies (e.g., Ref. 33) have shown that, in the absence of gravitational effects, there is very little difference in the average behavior of 2D and 3D simulations.

We shall analyze finger growth, the number of fingers observed in a displacement as a function of time and their average behavior in a first contact miscible displacement through a semi-infinite, two-dimensional (2D) rectilinear, porous medium (Fig. 1). The system is initially filled with a fluid of viscosity µ1. A lower viscosity fluid (viscosity µ2) that is first contact miscible with the ambient fluid is injected continuously, at constant rate, into the system along the left-hand boundary. There are no-normal flow boundaries top and bottom. The system is horizontal, so gravity can be neglected. It is assumed that both the fluids and the porous medium are incompressible.

FIG. 1.

The two-dimensional semi-infinite model studied. Fluid 2 is injected across the left hand face at constant rate. c is the mass fraction of the injected fluid. xD and yD are the dimensionless distances in x and y directions, respectively. There are no-normal flow boundaries at top and bottom.

FIG. 1.

The two-dimensional semi-infinite model studied. Fluid 2 is injected across the left hand face at constant rate. c is the mass fraction of the injected fluid. xD and yD are the dimensionless distances in x and y directions, respectively. There are no-normal flow boundaries at top and bottom.

Close modal

If the two fluids are first contact miscible, then the equation for the conservation of volume of the injected fluid is given by

ϕct+vc=ϕDc,
(1)

where c is the mass fraction of the injected fluid, ϕ is the porosity, D is the dispersion tensor, ∇ is the vector differential operator, and v is the Darcy velocity given by

v=KμmP,
(2)

where K is the permeability of the porous medium, µm(c) is the mixture viscosity, and P is the fluid pressure. Lohrenz et al.34 showed that the quarter power mixing rule is appropriate for miscible hydrocarbon mixtures in oil reservoirs

μm=1cμ11/4+cμ21/4.
(3)

Many authors analyzing viscous fingering have found it more convenient analytically to use the exponential mixing rule,

μm=μ1ec ln M,
(4)

where M=μ1μ2. Practically, these different functional forms give very similar mixture viscosities (Fig. 2).

FIG. 2.

Mixture viscosity as a function of concentration (µ1 = 10, µ2 = 1) obtained for two different functions used in the literature [Eqs. (3) and (4)].

FIG. 2.

Mixture viscosity as a function of concentration (µ1 = 10, µ2 = 1) obtained for two different functions used in the literature [Eqs. (3) and (4)].

Close modal

In general, the dispersion tensor, D, is velocity dependent (see, for example, the work of Perkins and Johnston50) but not dependent on concentration, and for linear flow, it is written as

D=DL00DT=Do+αLv00Do+αTv,
(5)

where Do is the molecular diffusivity in the porous medium, αL is the longitudinal dispersion coefficient, αT is the transverse dispersion coefficient, and v is the local speed. This velocity dependence of dispersion can complicate the analysis and interpretation of behaviors, particularly when we want to consider the impact of anisotropic dispersion, as it may result in more dispersion in the leading part of fingers (which are moving more quickly) and less in the trailing interfingered zones (which are moving more slowly). Zimmerman and Homsy21 simplified their analysis by assuming that only longitudinal dispersion was velocity dependent and that transverse dispersion only depended on the molecular diffusivity in the porous medium. In this paper, we neglect the velocity dependence of both components of dispersion although we consider the effects of different longitudinal and transverse diffusion.

Equations (1) and (2) are constrained by the equation for the conservation of total volume (as the system is assumed to be incompressible),

v=0.
(6)

Previous workers have suggested a number of different ways of scaling the system. In particular, we note that Tan and Homsy7 and Zimmerman and Homsy11 scale all lengths by the longitudinal dispersive length (DL/|v|) as their system has periodic boundaries and hence has no natural physical length scale. In contrast, Nijjer et al.22 choose the thickness of the domain as the reference physical length scale. As the system studied in this paper is semi-infinite and bounded perpendicular to flow, we follow Nijjer et al.22 in choosing the reference length scale to be the system width, H, so the dimensionless distances are

xD=xH,
(7)
yD=yH,
(8)

where x is the physical distance measured parallel to the average direction of flow and y is the physical distance measured transverse to the average direction of flow. We define the dimensionless time to be

tD=tDTH2
(9)

although we present some of our results in terms of the dimensionless advective time

tD,A=tvtHϕ.
(10)

The choice of tD as the dimensionless time is motivated by the analysis of Peters et al.35 and Tan and Homsy,7 which showed that transverse dispersion is the main control on initial finger number and growth.

Neglecting the velocity dependence of dispersion allows us to rewrite Eq. (1) as

ϕctD+PeTϵvDDc=D100ϵDc,
(11)

where ϵ = DT/DL, PeT=HvtϕDT is the transverse Peclet number, vD=vvt, and vt=1H0Hvxdy is the average speed in the direction of flow.

Although viscous fingering is highly nonlinear, in general, at early times, various authors (e.g., Refs. 1, 18, and 8) have shown that it is possible to use perturbation theory to analyze early time behavior. This involves writing the primary variables (c, v, P) in the form

a=āxD,tD+axD,yD,tD,
(12)

where ā is the mean value of the quantity (averaged transverse to the flow) and a′ is the fluctuation to the mean such that

01adyD=0
(13)

and is sufficiently small that all higher order terms can be neglected. Most authors assume that the perturbations take the form of a periodic fluctuation of a wavenumber k described by

a=ΦxD,tDeikyD,
(14)

where ΦxD,tD=ΦxD,0eωtD, i.e., the amplitude of the fluctuation grows at a rate of ω and the base state, ā, is a step function in concentration at tD = 0. They also assume periodic boundary conditions.

Tan and Homsy,8 amongst others, showed that there is a dimensionless finger wavenumber km with maximum dimensionless growth rate ωm and a cutoff wavenumber kc above which the perturbation will not grow. This cutoff wavenumber depends on diffusion/dispersion: if there were no diffusion/dispersion, then growth rate would grow linearly with the wavenumber and there would be no cutoff wavenumber (e.g., see analysis of Christie and Bond12). If the viscosity varies exponentially with concentration [Eq. (4)], then Tan and Homsy8 were able to obtain an analytical expression for the initial growth rate (at tD = 0) as a function of wavenumber, mobility ratio, and dispersion given by

ω=12PeT2ϵkRPeTϵkk2+2RPeTϵ+21ϵk,
(15)

where R = ln M and we have converted their equation into our system of dimensionless units.

Solving this equation for the cutoff wavenumber which occurs when ω = 0 gives

kc=R2PeTϵ1ϵ+ϵ.
(16)

The most dangerous wavenumber km occurs when dωdk=0. In this case, Eq. (15) can only be solved exactly when diffusion is isotropic, giving

km=254RPeT40.118RPeT
(17)

and

ωm=42+573535RPeT420.0225R2PeT2.
(18)

Tan and Homsy8 provide two approximate solutions for when ϵ ≫ 1 (transverse diffusion is much greater than longitudinal),

kmRPeT41322ϵ,
(19)
ωm=2ϵRPeT421212ϵ,
(20)

and when ϵ ≪ 1 (the more usual case, when transverse diffusion is much less than longitudinal),

km=2ϵRPeT4ϵ1353,
(21)
ωm=4RPeT4213ϵ13.
(22)

We see that in all cases, the wavenumber scales proportionally to (RPeT) and the growth rate of the fingers scales proportionally to (RPeT)2. Note that in our system of dimensionless parameters, the dimensionless wavenumber is k = 2πnf, where nf is the number of fingers. Manickam and Homsy36 showed that the preceding analysis can be extended to any function relating mixture viscosity to concentration by replacing R with Λ, where Λ is defined as

Λ=dμmdcc=0+dμmdcc=1μ1+μ2.
(23)

Overall, we see from this analysis that the initial fingering pattern (number of fingers and growth rate) is controlled by the transverse dispersion and the viscosity ratio (as originally suggested by Slobod and Thomas37).

If fingers grow from an initial concentration profile that is itself changing in time (due to longitudinal diffusion), then a similar behavior is observed, but kc, km, and ωm reduce with time. This is because the concentration gradient in the direction of flow is reducing due to longitudinal diffusion and the displacement becomes more stable over time. It is not possible to obtain an analytical solution for wavenumber and growth rates in this case although Tan and Homsy8 were able to obtain numerical solutions for the case R = 3 and ϵ = 1. Figure 3 sketches the overall form of Tan and Homsy’s8 results, while Fig. 3(b) shows the number of fingers corresponding to the most dangerous wavenumber (nf) as a function of dimensionless advective time, tD, A. These data were determined from Fig. 3 in the work of Tan and Homsy.8 The number of fingers decreases rapidly initially and quickly asymptotes to a constant, much lower number, as does the growth rate. Table I lists the number of fingers corresponding to the most dangerous wavenumber calculated from Eq. (17) for R = 3 together with those obtained for tD, A = 0 and tD, A = 0.04 from Fig. 3 in the work of Tan and Homsy.8 The close agreement between the number of fingers predicted by Eq. (17) and digitized from Tan and Homsy’s8 Fig. 3 demonstrates the accuracy of the digitization and conversion of the data into our system of dimensionless parameters. The growth rates predicted from Eq. (18) and read off Fig. 3 in the work of Tan and Homsy8 are also given. The number of fingers at tD, A = 0.04 is approximately half that at tD, A = 0, while the growth rate is reduced by a factor of 4.

FIG. 3.

The perturbation analysis of Tan and Homsy8 shows that the initial number of fingers decreases as the initial sharp interface between the injected and displaced fluid diffuses over time. (a) Schematic plot of initial growth rate of viscous fingers as a function of wavenumber, showing the growth rate in the absence of diffusion and how longitudinal diffusion reduces the growth rate of higher wavenumber fingers. This is based on quantitative plots provided in the work of Tan and Homsy.8 (b) The initial number of fingers with the maximum growth rate as a function of dimensionless time, for various PeT, ϵ = 1, R = 3. Data are taken directly from Fig. 3 in the work of Tan and Homsy.8 The symbols correspond to data read off from Fig. 3 in the work of Tan and Homsy,8 while the lines are to guide the eye.

FIG. 3.

The perturbation analysis of Tan and Homsy8 shows that the initial number of fingers decreases as the initial sharp interface between the injected and displaced fluid diffuses over time. (a) Schematic plot of initial growth rate of viscous fingers as a function of wavenumber, showing the growth rate in the absence of diffusion and how longitudinal diffusion reduces the growth rate of higher wavenumber fingers. This is based on quantitative plots provided in the work of Tan and Homsy.8 (b) The initial number of fingers with the maximum growth rate as a function of dimensionless time, for various PeT, ϵ = 1, R = 3. Data are taken directly from Fig. 3 in the work of Tan and Homsy.8 The symbols correspond to data read off from Fig. 3 in the work of Tan and Homsy,8 while the lines are to guide the eye.

Close modal
TABLE I.

The influence of longitudinal diffusion on the number of fingers and growth rate over advective time at tD,A = 0 and tD,A = 0.04 for different PeT. R = 3 and ϵ = 1 in all cases. The number of fingers and growth rate predicted by Eqs. (17) and (18) are shown together with the values read off Fig. 3 in the work of Tan and Homsy.8 The number of fingers at tD,A = 0.04 is approximately half that at tD,A = 0, and the growth is approximately one quarter that without diffusion included.

PeT, ϵ = 1nf, Eq. (17)nf, tD,A = 0nf, tD,A = 0.04nf,0/nf,0.04ωm Eq. (18)ωmtD,A = 0.04
100 5.6 5.6 3.7 1.5 22.5 6.6 
200 11.3 11.1 6.1 1.8 45 13.2 
300 16.9 16.7 7.8 2.1 67.5 19.7 
400 22.5 22.3 10.9 2.0 90 26.3 
500 28.2 27.9 12.7 2.2 112.5 32.9 
PeT, ϵ = 1nf, Eq. (17)nf, tD,A = 0nf, tD,A = 0.04nf,0/nf,0.04ωm Eq. (18)ωmtD,A = 0.04
100 5.6 5.6 3.7 1.5 22.5 6.6 
200 11.3 11.1 6.1 1.8 45 13.2 
300 16.9 16.7 7.8 2.1 67.5 19.7 
400 22.5 22.3 10.9 2.0 90 26.3 
500 28.2 27.9 12.7 2.2 112.5 32.9 

Let us now investigate the scaling of the growth rate and number of fingers at early time. Using the scaling described earlier, the number of fingers is related to the wavenumber km by

nf=km2π.
(24)

Examining Eqs. (17)(19), we see that in all cases,

kmRPeT.
(25)

Motivated by this observation, we divided both sides of Eq. (15) (the equation describing the growth rate of the most dangerous wavenumber at early time) by R2 to obtain

ωmR2=12PeT2ϵkmR1PeTϵkmRkmR2+21PeTϵkmR+21ϵkmR.
(26)

We can then find an expression for the scaled maximum wavenumber kmR by solving dωm/R2dkm/R=0 numerically. The results are shown in Fig. 4 where we plot km2πR=nfR as a function of PeT for ϵ = 1 × 10−5, 0.2, 1, 5, 1 × 105. We also plot the analytical solutions for ϵ = 1 [Eq. (17)] and ϵ ≫1 [Eq. (19)] both of which match the numerical solutions.

FIG. 4.

The maximum number of fingers at t = 0 predicted by the analytical solution obtained using linear stability analysis [Eq. (27)], assuming a sharp interface between displaced and displacing fluid (Ref. 8).

FIG. 4.

The maximum number of fingers at t = 0 predicted by the analytical solution obtained using linear stability analysis [Eq. (27)], assuming a sharp interface between displaced and displacing fluid (Ref. 8).

Close modal

For any given ϵ, there is a linear relationship between nfR and PeT. We also see that the greater ϵ is, the greater nfR is. However, there is an upper limit to this relationship because [from Eq. (19)] when ϵ → ∞,

nfR=18πPeT,
(27)

however, when ϵ → 0 (corresponding to a very large longitudinal diffusion), there are no fingers, as expected from Eq. (21), because the front is completely smeared out.

There are four fingering regimes that can be observed as miscible viscous fingers develop over time:

  1. Very early time, during which no fingers are observed because the growth rate of the fingers is less than the spreading (SP) of the interface due to longitudinal diffusion and dispersion. This seems to have been first described by Perkins et al.38 although it was subsequently also described by Kempers,49 Bacri et al.,39 and Loggia et al.40 

  2. Early time, when the growth rate of the fingers outstrips the rate of longitudinal mixing. The number and growth rate of the fingers can be calculated using linear perturbation theory such as that of Tan and Homsy,8 described in Sec. II D.

  3. Intermediate time when the fingers begin to interact nonlinearly. Depending on the viscosity ratio, the anisotropy in diffusion, and the Peclet number, the fingers may branch, merge, or fade although ultimately it seems that merging and fading dominate, so the total number of fingers reduces.

  4. Late time, when a single finger is left and propagates through the system. This was first described by Saffman and Taylor1 and analyzed by Outmans41 who showed that the finger width occupies approximately half of the overall system width. The transition and duration of this regime has been analyzed more recently by De Wit et al.10 and Nijjer et al.22 

The intermediate time regime and the transition to the late time regime are of particular interest as these are the regimes most often encountered in physical systems of engineering interest such as CO2 sequestration and enhanced oil recovery. Although the number of fingers may increase at specific time instances through splitting, it seems that on average, in miscible displacements in porous media, the number of fingers decreases with time. Splitting tends to occur at the finger tips for lower viscosity ratios (e.g., M = 47.542) and along the finger for higher viscosity ratios (e.g., M = 38343). It seems that the Peclet number in laboratory scale experiments is such that diffusion and dispersion tend to limit the ability of fingers to split. Perkins et al.38 described the overall reduction in the number of fingers with time as being due to the mechanisms of merging and fading, ascribing the fading of trailing fingers as being due to the growth of leading fingers which grow slightly more quickly and thus take more of the lower viscosity solvent. In contrast, Zimmerman and Homsy21 explained the fading of fingers as being due to spreading of the leading finger and thus blocking or adsorbing the trailing finger.

At late times, when only one single finger remains, there remains some uncertainty as to the ultimate fate of the finger. Tan and Homsy7 observed, in numerical simulations, that the single finger would split, one of the two fingers would fade forming a single finger and then this process would repeat. Zimmerman and Homsy11 found that this was less likely to occur with anisotropic dispersion (when the transverse dispersion was weak). They speculated that the behavior of this finger at late times would be independent of the Peclet number and the ratio of transverse to longitudinal dispersion for large Peclet numbers. Interestingly, Nijjer et al.22 do not see any tip splitting (TS) in the late time, single finger regime despite using isotropic diffusion. Furthermore, their analysis suggests that once the single finger forms, it decays exponentially due to transverse diffusion.

The dynamics of fingering during the intermediate and late time fingering regimes cannot be described analytically, unlike the very early and early time behaviors. They can be investigated using numerical simulation, but this investigation can be challenging due to the large system aspect ratios required and the fine meshes needed to minimize numerical diffusion.

In this paper, we find numerical solutions to Eqs. (1)–(5) using a finite difference, implicit-pressure explicit-saturation (IMPES) simulator developed by Christie and Bond44 to investigate the finger dynamics across all four regimes. No-normal flow boundary conditions are imposed along the boundaries parallel to the principal direction of flow. The simulator uses a quarter power mixing rule for viscosity. Second order accuracy is achieved using two-point, upstream weighting in combination with flux corrected transport (Ref. 45) to prevent unphysical oscillations around the shock front in immiscible displacements. The numerical calculations use double precision. The ability of the simulator to predict viscous fingering in 2D miscible displacements has previously been validated by comparison with pseudo 2D experimental results obtained from glass bead-packs by Christie and Bond,12 Christie et al.,46 and Davies et al.,47 amongst others.

A range of system aspect ratios were investigated from square to L/H = 10, together with a range of viscosity ratios, transverse Peclet numbers, and diffusion anisotropies (Table II). The investigations used anisotropic diffusion rather than velocity dependent dispersion to ensure that ϵ was constant across the model and thus simplifies analysis of the results. Further work is needed to compare these results with those that would be obtained if only longitudinal dispersion is velocity dependent (as assumed by Zimmerman and Homsy21) and if both transverse and longitudinal dispersion are velocity dependent, as is usually assumed to be the case on the field scale. As noted earlier in this paper, velocity dependent dispersion would lead to higher dispersion around a finger compared to the rear of the fingered front.

TABLE II.

Summary of data used in the investigation of fingering regimes.

Very early toIntermediate
Parametersintermediate timesto late time
Domain aspect ratio, (L/H10 
No. of grid blocks 600 × 200 2000 × 200 
Viscosity ratio, M 20, 55, 148 (R = 3, 4, 5) 20 (R = 3) 
Transverse Peclet 125, 250, 333, 500 125, 250, 500 
number, PeT   
Ratio of transverse to 0.2, 1, 5 0.25, 0.5, 1, 2 
longitudinal dispersion, ϵ   
Very early toIntermediate
Parametersintermediate timesto late time
Domain aspect ratio, (L/H10 
No. of grid blocks 600 × 200 2000 × 200 
Viscosity ratio, M 20, 55, 148 (R = 3, 4, 5) 20 (R = 3) 
Transverse Peclet 125, 250, 333, 500 125, 250, 500 
number, PeT   
Ratio of transverse to 0.2, 1, 5 0.25, 0.5, 1, 2 
longitudinal dispersion, ϵ   

The simulation domain was discretized into rectangular grid blocks, and the number of grid blocks was chosen using a grid refinement study to ensure that physical diffusion dominated over numerical diffusion. The simulator calculates each time step using a Courant-Friedrichs-Lewy (CFL) criterion to minimize numerical diffusion (Ref. 48). In this study, we chose CFL = 0.4 (slightly less than the theoretical value of 0.5 for a second order scheme) to allow for a small margin for error in the simulator estimates of velocity while minimizing numerical diffusion due to time stepping.

Miscible solvent was injected at a constant rate across the inlet face. This was matched by the rate of production from the outlet face. The fingers were explicitly triggered by using a uniformly distributed random concentration (0 < c < 1) at the inlet when t = 0, c = 0 elsewhere at initial time. A different, random initialization was used in every simulation.

We performed two sets of simulations, one to investigate the very early to intermediate time regime and the second to investigate the intermediate to late time regime. The first set used a square domain and a 600 × 200 grid so ∆xD = 1/600 = 0.001 67 and ∆yD = 1/200 = 0.005. The higher resolution was used in the x direction to ensure that the simulation captured the early time growth rates of the fingers more accurately. The second set of simulations required a longer domain; hence, we used L/H = 10 and a grid of 2000 × 200, giving ∆xD = ∆yD = 1/200 = 0.005. A comparison of simulations using these square and long aspect ratio domains showed that the different mesh resolution in the flow direction did not alter the fingering dynamics during the transition from early to intermediate times. The simulation input parameters used are summarized in Table II.

The simulation outputs were analyzed in terms of average concentration profiles (averaged transverse to the principle flow direction), growth of fingers, and the number of fingers vs dimensionless advective time, tD,A (pore volumes injected for a unit aspect ratio system), given by

tD,A=tvtHϕ,
(28)

as well as dimensionless time based on transverse diffusion [Eq. (9)].

To quantify the growth of the fingers, we express their spreading in the direction of flow in terms of the dimensionless mixing length Lmix,

Lmix=xDc=0.02xDc=0.98,
(29)

where xD|c=0.02 and xD|c=0.98 are the dimensionless distances in the direction of flow where the average concentration equals c = 0.02 and c = 0.98, respectively.

We used contour lines to keep track of the evolution of the fingers as a function of time. This was necessary because of the finite difference simulation approach adopted. In spectral methods, this can be done by directly taking the wavenumber with highest power, e.g., Zimmerman and Homsy.21 At every 5th time step, we computed the contour of the concentration at a specified concentration level clevel using the marching square method. Due to coalescence of the fingers and blob formation due to pinchoff of the fingers, there could be multiple unconnected contour lines for a given clevel as illustrated in Fig. 5(a). In this instance, only the longest contour line would be selected (and the rest of the lines ignored). We then found the local minima and maxima along this selected contour line which correspond to the peaks and troughs. We excluded any peak if it was located less than two grid blocks further toward the producer than its adjacent trough, as shown in Fig. 5(b). This was particularly important at very early times to exclude the initial perturbation introduced at t = 0. clevel = 0.5 was used in this study. This value was found to capture the number of fingers correctly, as shown in Sec. V, but sensitivity studies showed that results did not change significantly if a different contour was analyzed.

FIG. 5.

Illustration of the number of fingers nf was calculated using the contour lines for c = 0.5. (a) Only the longest contour line was selected for calculation of nf; (b) nf is given by the number of peaks. We exclude any peak that is located less than two grid blocks further toward the producer than its adjacent trough.

FIG. 5.

Illustration of the number of fingers nf was calculated using the contour lines for c = 0.5. (a) Only the longest contour line was selected for calculation of nf; (b) nf is given by the number of peaks. We exclude any peak that is located less than two grid blocks further toward the producer than its adjacent trough.

Close modal

We also tracked the trajectory of each finger peak in time in order to investigate the dynamics of the fingers during the intermediate to late time regimes.

Before investigating different aspects of finger growth and decay, we provide a description of the overall finger dynamics using the different methods applied to analyze them. Figure 6 shows the growth and decay of fingers predicted by numerical simulation for PeT = 250, ϵ = 1, and R = 3. Figure 6(a) shows the number of fingers nf as a function of tD,A, Fig. 6(b) shows the corresponding mixing length against time on a log scale, and Fig. 6(c) provides snapshots of the corresponding concentration distributions at different times.

FIG. 6.

Illustration of the different methods used to analyze the fingering behavior using the method in Sec. IV B. (a) Number of fingers as a function of dimensionless advective time tA,D. There are four regimes observed here: I: very early time, II: early time, III: intermediate time, and IV: late time. (b) Mixing length of fingered zone against dimensionless advective time (both scales logarithmic) showing a diffusive growth at very early time before entering a linear growth at intermediate to late time regimes. (c) Concentration maps for each of the regimes (at tA,D = 0.05, 0.3, 1, 5) showing the typical displacement patterns associated with each regime.

FIG. 6.

Illustration of the different methods used to analyze the fingering behavior using the method in Sec. IV B. (a) Number of fingers as a function of dimensionless advective time tA,D. There are four regimes observed here: I: very early time, II: early time, III: intermediate time, and IV: late time. (b) Mixing length of fingered zone against dimensionless advective time (both scales logarithmic) showing a diffusive growth at very early time before entering a linear growth at intermediate to late time regimes. (c) Concentration maps for each of the regimes (at tA,D = 0.05, 0.3, 1, 5) showing the typical displacement patterns associated with each regime.

Close modal

The four different fingering regimes can be clearly seen in Fig. 6(a). At very early time, there are no fingers and the mixing length of the interface grows diffusively [Fig. 6(b)] until at around tD,A = 0.3, 8 small fingers form. The number of fingers then decays monotonically until at around tD,A = 0.9, only 3 fingers remain. Just after tD,A = 1, a tip splitting event occurs [from Fig. 6(c), we see that this occurs in the leading finger]. The number of fingers then slowly decreases, apart from at around tD,A = 2.3 when another splitting event occurs, until just before tD,A = 5 only one finger remains, and the system enters the late time regime. Figure 6(b) shows that the mixing length of the fingered zone grows linearly with time throughout the intermediate time regime. This is consistent with results obtained by many previous researchers (e.g., Refs. 28 and 29) who have shown that fractional flow based empirical models such as those of Koval25 and Todd and Longstaff26 can be used to estimate the average behavior of unstable viscous fingering in linear systems.

Having demonstrated that we can predict and analyze the full range of fingering regimes in our simulations, we now go on to investigate the different fingering regimes in more detail together with the scaling of the transitions between those regimes.

As noted in Sec. III, at very early time, the sharp front between the displacing and displaced fluids is diffused by longitudinal dispersion DL. This initially spreads the front in the direction of flow more quickly than any fingers can grow due to the very high concentration gradient across the front.

Figure 7(a) shows the influence of R and PeT on the transition between diffusive flow at very early time to viscous-dominated unstable flow, in terms of mixing length vs advective time, when ϵ = 1. For a given PeT, Lmix at very early time does not depend upon R; however, the crossover from diffusion to fingering is earlier for higher R. This is as expected from linear perturbation analysis as the early time growth rate of the fingers increases approximately as R2 [Eq. (30)].

FIG. 7.

The crossover between the very early time to early time regime. (a) Mixing lengths for R = 3, 4 and PeT = 250, 500, ϵ = 1 calculated from detailed simulations. (b) The crossover times estimated for R = 3, PeT = 500, ϵ = 1 by comparing the longitudinal spreading due to DL at very early time [Eq. (30)] with: (i) the Koval model (as suggested by Kempers49) (ii) the growth rate predicted by Eq. (8) (as suggested by Bacri et al.39), and (iii) the growth rate predicted by Tan and Homsy8 for tD,A = 0.04 (here, we used the initial longitudinal disturbance length lD of 1 grid block in the direction of flow).

FIG. 7.

The crossover between the very early time to early time regime. (a) Mixing lengths for R = 3, 4 and PeT = 250, 500, ϵ = 1 calculated from detailed simulations. (b) The crossover times estimated for R = 3, PeT = 500, ϵ = 1 by comparing the longitudinal spreading due to DL at very early time [Eq. (30)] with: (i) the Koval model (as suggested by Kempers49) (ii) the growth rate predicted by Eq. (8) (as suggested by Bacri et al.39), and (iii) the growth rate predicted by Tan and Homsy8 for tD,A = 0.04 (here, we used the initial longitudinal disturbance length lD of 1 grid block in the direction of flow).

Close modal

Kempers49 argued heuristically that the time for this crossover can be roughly estimated by equating the mixing length induced by longitudinal diffusion,

Lmix,D=5.8tD,APeL,
(30)

to the linear spreading based on the Koval25 model,

Lmix,K=Me1MetD,A,
(31)

where Me is the effective viscosity ratio defined as

Me=0.78+0.22M144.
(32)

Note that we have adapted the Kempers49 equation for diffusion [Eq. (30)] so that it calculates the mixing length from c = 0.02 to c = 0.98, so it is consistent with the mixing length defined in Eq. (29), while the mixing length predicted by the Koval25 model, Eq. (31), is from c = 0 to c = 1.

Bacri et al.39 proposed that this crossover can be estimated by equating the mixing length from longitudinal diffusion [Eq. (30)] to the growth of the initial perturbation (without diffusion) estimated from linear stability theory,

Lmix,P=lDeωmtD,
(33)

where lD is the initial perturbation length in the x direction.

We plot the growth of the mixing length using Eqs. (31)(33) together with the prediction from detailed numerical simulation in Fig. 7(b) for R = 3, PeT = 500, and ϵ = 1. Here, lD was taken as the dimensionless length of 1 grid block in the x direction as this is the length of the initial random perturbation made in the simulation. We also show the effect of using the growth rate predicted by Tan and Homsy8 at tD,A = 0.04 in Eq. (33), when the front is smeared by longitudinal diffusion. The intersection of the diffusive mixing length (plotted with a solid black line) with either the mixing length vs time calculated from Eq. (31), or Eq. (33), indicates the predicted crossover from the very early time to the early time regime for each model. It seems that both the approaches predict a much earlier transition to viscous fingering than is observed in the simulation. The Kempers49 approach predicts a very early transition probably because the Koval model was tuned to estimate average finger growth in the late time regime when the mixing length grows linearly. The Bacri et al.39 approach using the analytical calculation of the growth rate predicts a slightly later transition (at tD,A ∼ 0.04) than the Kempers49 approach, but this is considerably earlier than the transition seen in the simulation (which transitions to early time fingering at around tD,A = 0.15). This is probably because longitudinal diffusion reduces the growth rate of the fingers, as using the growth rate predicted by Tan and Homsy8 for tD,A = 0.04 results in a much later estimate of a transition to viscous fingering (at tD,A ∼ 0.09). Even this growth is probably at overestimate of the actual growth rate as Table I shows that the growth rate will continue to reduce at larger values of tD,A.

Figure 8 shows the number of fingers obtained at early time from simulation, scaled by R, as a function of PeT and ϵ. These follow the same linear trend as the maximum number of fingers predicted from the perturbation analysis of Tan and Homsy8 as shown in Fig. 4. The number of fingers increases with the ratio of transverse to longitudinal dispersion for a given PeT because there is less smoothing of the initial step concentration profile by longitudinal dispersion at early time. There are fewer fingers for lower ϵ because there is more smearing of the front by longitudinal dispersion at early time and hence the front is less unstable.

FIG. 8.

Maximum number of fingers obtained from the simulation (see Fig. 10) for R = 3, 4, 5 compared to the analytical solutions for the maximum number of fingers predicted by linear stability analysis [Ref. 8, Eq. (13)] assuming an exponential mixing rule for viscosity. The maximum number of fingers obtained in the simulations is approximately half that predicted analytically.

FIG. 8.

Maximum number of fingers obtained from the simulation (see Fig. 10) for R = 3, 4, 5 compared to the analytical solutions for the maximum number of fingers predicted by linear stability analysis [Ref. 8, Eq. (13)] assuming an exponential mixing rule for viscosity. The maximum number of fingers obtained in the simulations is approximately half that predicted analytically.

Close modal

It can be seen, however, that the gradients of nf vs PeT in Fig. 8 are only half of those predicted in Fig. 4, i.e., our simulations only give half the number of fingers predicted by the analytical solutions from linear stability theory. This reduction is consistent with the analysis of Tan and Homsy8 when the effect of longitudinal diffusion smearing out the front is taken into account (see Table I). This is not due to the simulations using a quarter power mixing rule, rather than an exponential mixing rule. Figure 9 shows the theoretical number of fingers expected using Eq. (23) in Eq. (17) (which neglects the effect of longitudinal diffusion on the number of fingers) for the quarter power mixing rule for viscosity (which is used in the simulations), the exponential mixing rule and the linear mixing rule. This shows that using quarter power mixing for viscosity should result in even more fingers than using an exponential mixing rule.

FIG. 9.

The influence of the viscosity mixing function on the maximum number of fingers predicted by the linear stability analysis of Tan and Homsy8 for isotropic diffusion when M = 148.

FIG. 9.

The influence of the viscosity mixing function on the maximum number of fingers predicted by the linear stability analysis of Tan and Homsy8 for isotropic diffusion when M = 148.

Close modal

This reduced number of initial fingers is also seen in other studies. Suekane et al.42 found that the maximum number of fingers measured in their laboratory experiments was approximately two orders of magnitude fewer than the number predicted by linear stability analysis. Zimmerman and Homsy21 too found that the number of fingers formed initially in their simulations was often half that predicted by linear stability theory. They used a simulator with periodic boundary conditions. They explained the reduction in the number of fingers by a pairwise interaction mechanism rather than slower growth and a reduced number of fingers due to front smearing. We cannot find any other comparisons of numbers of fingers predicted by linear stability analysis with those obtained either by numerical simulation or in the laboratory. We note that Tchelepi et al.51 compared the growth rates observed experimentally with those predicted by linear stability analysis and found that the experimental growth rates were much larger than those predicted analytically. They suggested that this was because the finger dynamics in their experiments were nonlinear even at the early times observed.

We suggest that the reduced number of fingers/reduced growth rate observed in all our simulations at very early time is consistent with longitudinal diffusion smearing out the initial step change in concentration between displaced and displacing fluid and reducing the instability of the displacement, as discussed by Heller52 and analyzed by Tan and Homsy.8 This hypothesis should be further verified by performing a series of pseudo 2D beadpack experiments in which the initial number of fingers and their early growth are monitored.

We now discuss some of the ways that fingers interact in the intermediate time regime with reference to the results shown in Fig. 10, where R = 3, 4, 5, and PeT = 250, 500, and ϵ = 1. The figure shows the concentration distribution at tD,A = 0.5 alongside the evolution of the 50% concentration contour up to that time. The finger tips are marked with different colors to indicate the dimensionless advective time tD,A. As expected, initially more fingers start growing sooner and grow more quickly for larger PeT or R. The fingering dynamics become nonlinear sooner with more splitting/branching occurring at higher viscosity ratios and PeT.

FIG. 10.

The early to intermediate time dynamics of the fingers (tD,A = 0.5) with isotropic dispersion for PeT = 250, 500 and R = 3, 4, 5 (M = 20, 55, 148). The number of fingers and the interaction between the fingers increase as we increase R and PeT. We observe various events including spreading (SP), shielding (SH), coalescence (CO), double coalescence (DC), tip splitting (TS), and side branching (SB).

FIG. 10.

The early to intermediate time dynamics of the fingers (tD,A = 0.5) with isotropic dispersion for PeT = 250, 500 and R = 3, 4, 5 (M = 20, 55, 148). The number of fingers and the interaction between the fingers increase as we increase R and PeT. We observe various events including spreading (SP), shielding (SH), coalescence (CO), double coalescence (DC), tip splitting (TS), and side branching (SB).

Close modal

Examining Fig. 10, we can identify a number of nonlinear fingering behaviors including

  1. Shielding (abbreviated as SH in Fig. 10), which occurs when a finger that is slightly ahead of its neighbors becomes more dominant and consequently halts the growth of its neighboring fingers.

  2. Spreading (SP), which is the lateral growth of the fingers which reduces the average wavelength of the fingers. This can either be purely due to transverse dispersion or as a result of finger coalescence. In Fig. 10, the spreading appears to be purely due to transverse diffusion.

  3. Coalescence (CO), which occurs when the fingers merge to become a single larger finger. In Fig. 10, this seems to result when a finger makes a sharp turn to coalesce with its neighbor. There are also instances of double coalescence (DC) as discussed by Islam and Azaiez15 where a pair of fingers adjacent to a dominant, shielding finger bend and merge into its base.

  4. Tip splitting (TS), which refers to the instability when a finger-tip splits into two. Tan and Homsy7 argue that, for this to happen, the finger-tip needs to spread so that it is wider than two fingers of the most dangerous wavelength for the PeT and ϵ of the system and the longitudinal concentration gradient at the tip.

  5. Side branching (SB) is also observed particularly in the highest mobility ratio, highest PeT case.

In general, we observe that the frequency of splitting events increases when we have a higher mobility ratio and PeT.

Figure 11 shows the effects of increasing and reducing longitudinal dispersion (ϵ = 0.2, 5) when PeT = 250. As expected from Fig. 4, more fingers form and grow more rapidly for ϵ = 5 compared to the isotopic case, while fewer fingers form that grow more slowly for ϵ = 0.2. For a given PeT, increasing longitudinal dispersion stabilizes the flow with fewer fingers formed as shown for ϵ = 0.2, while decreasing it results in more fingers, as shown for the case when ϵ = 5. There are also more side branching events in this case as the flow is influenced more heavily by DT. Although this case is unrealistic (as longitudinal dispersion DL is typically higher than DT for flow in porous media), it provides further evidence that a smaller DL will result in a greater flow instability which is consistent with the early time linear stability analysis.

FIG. 11.

The effects of anisotropic dispersion (ϵ = 0.2, 5) on the early time growth of the fingers (tD,A = 0.5) for PeT = 250 and R = 3, 4, 5. The corresponding cases with isotropic dispersion are shown in Fig. 10.

FIG. 11.

The effects of anisotropic dispersion (ϵ = 0.2, 5) on the early time growth of the fingers (tD,A = 0.5) for PeT = 250 and R = 3, 4, 5. The corresponding cases with isotropic dispersion are shown in Fig. 10.

Close modal

Figure 12 shows the number of fingers, nf, as a function of tD,A for PeT = 500, ϵ = 1, and R = 3, 4, 5 (the corresponding concentration maps are shown in Fig. 10). The onset of fingering occurs earlier for R = 5 as discussed in Sec. V A. After the initial growth of the fingers, nonlinear mechanisms such as shielding and coalescence reduce the number of fingers relatively quickly. The laboratory experiments by Wooding53 showed that the number of fingers reduces diffusively as a function of time, i.e., nft−0.5, and this depended upon DT. Although we observed similar behavior for R = 3 for 0.1 < tD,A < 0.3 [see Fig. 11(b)], for higher mobility ratios, we observed the following:

  1. Initially, the number of fingers decreases approximately as nftD,A1 [Fig. 11(b)]. Nijjer et al.22 also observed this in their simulations and argued that there is a small period when advection dominates shortly after the early time regime. They proposed that nfRPeTtD,A1.

  2. After this initial period, for R = 4, 5, nf starts to increase as a result of branching and tip splitting, before there is further merging and coalescence. No tip splitting was observed in the R = 3 and PeT = 500 case for tD,A up to 0.5. This is similar to the observations by Zimmerman and Homsy21 who found that they had to increase PeT > 300 before tip splitting occurred in their spectral method simulations.

FIG. 12.

(a) The number of fingers from early to intermediate time for PeT = 500 and ϵ = 1. The corresponding concentration maps are shown in Fig. 10. (b) Same as (a) using logarithmic scales.

FIG. 12.

(a) The number of fingers from early to intermediate time for PeT = 500 and ϵ = 1. The corresponding concentration maps are shown in Fig. 10. (b) Same as (a) using logarithmic scales.

Close modal

Based on the results from early time perturbation analyses [Eqs. (16) and (17)] and dimensionless scaling described in Sec. II, we propose the following scaling for the number of fingers and time

nfRPeT,tD1RPeT2.
(34)

Note that tD is the dimensionless diffusive time defined in Eq. (9). Nijjer et al.22 proposed the same scaling for nf but scaled advective time by PeTR2 although they only considered isotropic diffusion and obtained the scaling based on physical arguments. They also used a logarithmic time in their graphs.

The results for different simulations and different scalings using PeT = 500, ϵ = 1, 0.2, 5, and R = 3, 4, 5 are shown in Fig. 13. We see that using the scaling in Eq. (34) [Fig. 13(b)] results in the curves collapsing into groups according to their ϵ values. This is because the scaling variables [Eqs. (17) and (18)] did not take into account the anisotropy of dispersion. Interestingly, scaling the number of fingers by the most dangerous wavenumber and diffusive time by the growth rate of that wavenumber [determined from Eq. (26)] results in all curves being superimposed suggesting that the nonlinear dynamics depend directly on the initially linear behavior of the fingers at early time [Fig. 13(c)]. This is despite the observation earlier that this linear analysis does not even predict the initial number of fingers and their growth rate correctly as it neglects the smearing of the front by longitudinal diffusion.

FIG. 13.

Different ways of scaling the number of fingers as a function of advective dimensionless time from very early to intermediate time regimes for PeT = 500, R = 3, 4, 5, and ϵ = 0.2, 1, 5. (a) Unscaled nf as a function of tD,A. (b) Scaling of nf using PeTR and tD,A by multiplying by R2PeT2 assuming isotropic diffusion [Eqs. (17) and (18)]. (c) Scaling of nf using the early the most dangerous wavenumber and growth rate determined from Eq. (26), which includes anisotropic diffusion.

FIG. 13.

Different ways of scaling the number of fingers as a function of advective dimensionless time from very early to intermediate time regimes for PeT = 500, R = 3, 4, 5, and ϵ = 0.2, 1, 5. (a) Unscaled nf as a function of tD,A. (b) Scaling of nf using PeTR and tD,A by multiplying by R2PeT2 assuming isotropic diffusion [Eqs. (17) and (18)]. (c) Scaling of nf using the early the most dangerous wavenumber and growth rate determined from Eq. (26), which includes anisotropic diffusion.

Close modal

We now investigate the dynamics and scaling of the intermediate to the late time fingering regime when the number of fingers reduces to one. This is the regime that is most likely to be observed in flows in the subsurface such as CO2 sequestration and enhanced oil recovery.

Figure 14 shows a typical set of results for the case R = 3, PeT = 500, and ϵ = 1. The aspect ratio of the simulation model is 10. Figure 14(a) shows how the 50% concentration contour evolves over dimensionless advective time. The corresponding concentration distribution shown in Fig. 14(b) shows the concentration map at tD,A = 5 for visual reference. Figures 14(c) and 14(d) show the number of fingers and mixing length as a function of tD,A, respectively. Note that the time scale in Fig. 14(d) is logarithmic.

FIG. 14.

Alternative ways of analyzing the results for R = 3, PeT = 500, and ϵ = 1. (a) The tracking of the finger tips using the method presented in Sec. IV B. The fingers tips are marked with different colors to represent advective time tD,A. (b) Concentration map at tD,A = 5. (c) Number of fingers as a function of advective time. (d) Mixing length as a function of advective time.

FIG. 14.

Alternative ways of analyzing the results for R = 3, PeT = 500, and ϵ = 1. (a) The tracking of the finger tips using the method presented in Sec. IV B. The fingers tips are marked with different colors to represent advective time tD,A. (b) Concentration map at tD,A = 5. (c) Number of fingers as a function of advective time. (d) Mixing length as a function of advective time.

Close modal

Figure 14(a) shows that the displacement starts with 12 fingers but very quickly (by tD,A = 1.8) nine of those fingers have faded due to shielding and coalescence. Two of the remaining fingers experience tip splitting at approximately tD,A = 2 and tD,A = 2.5, but there are only 2 fingers by tD,A = 3. In contradiction to the analysis of Yang and Yortsos24 we do not see any fingers that grow preferentially and more strongly along the boundaries. As seen in previous simulations, the mixing length grows proportionally to tD,A until fingers begin to grow at tD,A ∼ 0.1. After this time, the number of fingers decays monotonically (except for a brief period when there are tip splitting events between tD,A = 2 and tD,A = 3) to two fingers at tD,A = 3.3. The mixing length grows linearly with tD,A over this period.

We examine the scaling of the number of fingers and the mixing length during the transition to late time behavior in Fig. 15 for R = 3, PeT = 500, 250, 125, and ϵ = 0.25, 0.5, 1, 2. As before, we are able to collapse the plots for number of fingers against time by plotting nf/(km/2π) against tD,Aωm, as shown in Fig. 15(a).

FIG. 15.

The scaling of the number of fingers and mixing length during the intermediate to late time fingering regimes for R = 4, PeT = 125, 250, 500, and ϵ = 0.25, 0.5, 1, 2. (a) Number of fingers (divided by the most dangerous wavenumber) as a function of dimensionless advective time (multiplied by the growth rate of the most dangerous wavenumber), further confirming this scaling collapse all results into one curve; mixing length as a function of unscaled dimensionless time (b) isotropic diffusion; (c) PeT = 125, 250, 500 and ϵ = 0.25, 0.5, 1, 2; (d) PeT = 500, varying ϵ. In all cases, the early time mixing length scales as the intermediate to late time mixing length scales linearly with time.

FIG. 15.

The scaling of the number of fingers and mixing length during the intermediate to late time fingering regimes for R = 4, PeT = 125, 250, 500, and ϵ = 0.25, 0.5, 1, 2. (a) Number of fingers (divided by the most dangerous wavenumber) as a function of dimensionless advective time (multiplied by the growth rate of the most dangerous wavenumber), further confirming this scaling collapse all results into one curve; mixing length as a function of unscaled dimensionless time (b) isotropic diffusion; (c) PeT = 125, 250, 500 and ϵ = 0.25, 0.5, 1, 2; (d) PeT = 500, varying ϵ. In all cases, the early time mixing length scales as the intermediate to late time mixing length scales linearly with time.

Close modal

For isotropic dispersion, the mixing lengths Lmix for PeT = 500, 250, 125 exhibit different behaviors in the different fingering regimes, as shown in Fig. 15(b). A lower PeT gives a higher mixing length at very early time (as longitudinal dispersion DL is higher) but yields a lower mixing length at intermediate to late times.

Plotting Lmix against dimensionless advective time for simulations with the same longitudinal dispersion (DL = 0.004) but different PeT = 125, 250, 500 in Fig. 15(c) gives identical Lmix at very early time, as expected, however, at intermediate to late time regimes, Lmix depends mainly upon PeT, as shown in Fig. 15(d) where we plot the cases with the constant transverse Peclet number (PeT = 500) and different ϵ = 0.25, 0.5, 1.

Like Tan and Homsy,7 we see the tip splitting of the single remaining finger at intermediate to late times. Figure 16 compares the dynamics of the tip-splitting when PeT = 250 but with different longitudinal diffusivities (ϵ = 1, 2) as well as showing a case when tip-splitting did not occur (ϵ = 0.5). It also shows the mean concentration transverse to the principle flow direction as a function of time for each of these cases. When transverse diffusion is high (ϵ = 2), then the tip splitting process is repeated several times. The finger-tip splits and then one of the two fingers becomes dominant and the other fades away. The remaining dominant finger then splits again at its tip followed by one of the fingers fading.

FIG. 16.

The average concentration profiles for 0 < tD,A < 5 and the maps of finger dynamics at tD,A = 5. In all cases, PeT = 250. (a) ϵ = 2; (b) ϵ = 1; and (c) ϵ = 0.5. Reducing ϵ reduces the tendency of the flow to exhibit tip-splitting in the intermediate to late time regime.

FIG. 16.

The average concentration profiles for 0 < tD,A < 5 and the maps of finger dynamics at tD,A = 5. In all cases, PeT = 250. (a) ϵ = 2; (b) ϵ = 1; and (c) ϵ = 0.5. Reducing ϵ reduces the tendency of the flow to exhibit tip-splitting in the intermediate to late time regime.

Close modal

For the intermediate value of transverse diffusivity (ϵ = 1), only one such splitting event occurs while for the lowest value of transverse diffusivity no splitting event is seen. This is consistent with Tan and Homsy’s suggestion7 that the finger can only split when it is wide enough to sustain two fingers according to linear stability analysis as well as the observation of Zimmerman and Homsy21 using a simulator with periodic boundary conditions, that tip splitting was less likely to occur when ϵ < 1. As we have discussed previously, the most dangerous wavenumber depends upon the mobility ratio between the finger, the concentration gradient between the finger and the displaced fluid, and the displaced fluid and PeT. In the late time regime, the finger is formed of a mixture of displacing and displaced fluid (as shown in both the concentration distributions and graphs of mean concentration against tD,A in Fig. 16), and so there is a lower viscosity ratio between the finger and the displacing fluid than is the case initially. This is shown more clearly in Fig. 17 for the ϵ = 1 case, especially in Fig. 17(c) where we plot the maximum concentration in the finger against time. This means that the maximum number of fingers predicted by linear stability theory will be lower at later times, i.e., they will have a greater wavelength than at the beginning of the displacement. The maximum concentration in the finger-tip is controlled by the interplay of transverse dispersion, longitudinal dispersion, and viscosity ratio. A smaller DL results in a steeper concentration gradient in the direction of flow at the front of the finger, but a higher DT will tend to spread the finger transverse to flow. There is thus an optimum combination of DL, DT, and mobility ratio for tip-splitting. In particular, if DL is too large then there will be no tip-splitting [as observed in Fig. 16(c)]. This also means that tip-splitting is less likely to occur if ϵ < 1 as was also found by Zimmerman and Homsy.21 

FIG. 17.

Change in concentration along the path of the dominant finger in the late time regime over time for R = 3, PeT = 250, and ϵ = 1. (a) Illustration of the maximum concentration cmax,x (shaded black) along the dominant finger; (b) cmax,x (red lines) along the x axis over time (0 < tD,A < 5). We track the front of the finger-tip (marked as blue circles) by finding the concentration corresponding to a constant dcmax,x/dxD (chosen to be −10). (c) The concentration in the finger tip from (b) is gradually reducing suggesting there will be no tip-splitting at late time as the instability at the finger-tip will also be reducing.

FIG. 17.

Change in concentration along the path of the dominant finger in the late time regime over time for R = 3, PeT = 250, and ϵ = 1. (a) Illustration of the maximum concentration cmax,x (shaded black) along the dominant finger; (b) cmax,x (red lines) along the x axis over time (0 < tD,A < 5). We track the front of the finger-tip (marked as blue circles) by finding the concentration corresponding to a constant dcmax,x/dxD (chosen to be −10). (c) The concentration in the finger tip from (b) is gradually reducing suggesting there will be no tip-splitting at late time as the instability at the finger-tip will also be reducing.

Close modal

Tan and Homsy7 speculated that repeated tip-splitting could prevent the formation of a single dominant finger in some cases, whereas Nijjer et al.22 suggested that the system will always tend to form a single finger: they did not see any splitting of the single finger in their simulations (which used periodic boundary conditions), once formed. Our simulations and the discussion in the above paragraph suggest that tip-splitting, if it occurs, will delay the time at which a single finger forms but not prevent this. Figure 16 shows the average concentration along the path of the dominant finger over time for different values of ϵ. In Fig. 16(a) (where there are several tip-splitting events), the average concentration at the finger-tip is gradually reducing, suggesting that it will eventually decay to the point where the viscosity ratio between the tip and the displacing fluid is too low for tip-splitting to occur. This is seen for the cases shown in Figs. 16(b) and 16(c). We further note that the maximum concentration at the finger reduces according to tD,A once the finger-tip is no longer sufficiently unstable to split [Fig. 17(c)]. This suggests that the finger development is now controlled by transverse diffusion spreading the finger.

Figure 18 shows the same data from Fig. 15(a) but using a logarithmic scale for both axes. It confirms that once tip-splitting no longer occurs, the finger grows diffusively rather than advectively. We can use this scaling to estimate the beginning of the late time regime (when the remaining finger spreads diffusively and can no longer split). We notice from Fig. 18(b) that the decay of the final 2–3 fingers to a single dominant finger is described by nf2πkm1tDωm; thus, the time at km which all the fingers have disappear except for 1 can be estimated from

tD,11ωmkm2π2

(and there is only one finger left so nf = 1). Assuming isotropic diffusion and substituting for km from Eq. (17) and ωm from Eq. (18), we obtain

tD,10.016H2DT.
(35)
FIG. 18.

The same data as shown in Fig. 15(a) (the rescaled number of fingers against rescaled dimensionless diffusive time, for R = 3, PeT = 500, 250, 125, and ϵ = 0.25, 0.5, 2), but using a logarithmic scale on both axes. The rescaled number of fingers reduces diffusively (nf/km/2πtD,Aωm0.5) at intermediate to late times.

FIG. 18.

The same data as shown in Fig. 15(a) (the rescaled number of fingers against rescaled dimensionless diffusive time, for R = 3, PeT = 500, 250, 125, and ϵ = 0.25, 0.5, 2), but using a logarithmic scale on both axes. The rescaled number of fingers reduces diffusively (nf/km/2πtD,Aωm0.5) at intermediate to late times.

Close modal

The scaling of the number of fingers in the intermediate regime with 1/tD,A is consistent with the physical arguments presented by Nijjer et al.22 for isotropic diffusion although they did not discuss how to estimate the time at which the number of fingers decays to 1 or the impact of anisotropic diffusion.

Recasting Eq. (35) to obtain the dimensionless advective time, tD,A, at which the system decays to 1 finger results in

tD,A0.016HLPeT.
(36)

It is interesting that the time does not seem to depend directly upon R. Gardner and Ypma20 plotted unrecovered oil at tD,A = 1.1 from a series of physical experiments (performed in cores with L/H up to 20) and numerical simulations against the reciprocal of the dimensionless number,

NT=HLPeT.
(37)

The simulations and experiments were performed in different aspect ratio systems, yet they found that recovery scaled very well with 1/NT (see their Fig. 11). They found that the volume of unrecovered oil was constant for 1/NT < 0.02 and reduced as NT increased for 1/NT > 0.02. They suggested that this was because transverse diffusion stabilized the displacement; however, our analysis suggests that in fact the lower recovery occurs because there is only one finger in the system.

The lack of tip splitting seen at late time in the Nijjer et al.22 simulations could be due to the high Peclet numbers (1000 and 2000) used combined with isotropic diffusion. The above analysis suggests that it would take longer for a single finger to form in their simulations because of the high Peclet number. This finger would have a lower concentration at the tip because of the time taken to form and would also spread more slowly because of the low transverse diffusion. Consequently, for their values of R and PeT, the finger-tip may never be able to spread to the point that it can support two fingers before its concentration gradient is so low it is effectively stable. This argument suggests that tip splitting of the single finger may only occur for a small range of R and PeT; however, further simulations, supported by physical experiments, would be needed to substantiate this hypothesis.

This study has investigated the dynamics and scaling of miscible viscous fingering in a system with no flow boundaries perpendicular to the principle flow direction. This was performed using a second order finite difference simulator that has previously been validated against experiments. We have studied the influence of the transverse Peclet number, viscosity ratio, and the diffusion anisotropy on the behavior.

Overall the behavior is similar to that observed in previous studies (Refs. 21 and 22) that investigated fingering in systems with periodic boundary conditions. The main difference compared to previous studies was with Nijjer et al.22 We see tip-splitting after a single finger has formed, which they did not. We speculate that this is a result of isotropic diffusion and the high Peclet numbers (1000, 2000) in their late time simulations rather than the periodic boundary conditions as tip splitting was seen by Tan and Homsy8 and Zimmerman and Homsy.21 The combination of isotropic diffusion and high Peclet numbers means the finger-tip, for Nijjer et al.’s22 combination of R and PeT, never becomes sufficiently unstable to split. Our study investigated anisotropic diffusion and slightly lower Peclet numbers (PeT ≤ 500). Nonetheless, we observed that this tip splitting eventually stopped and, as suggested by Nijjer et al.,22 a single finger propagated through the displaced fluid. We did not see any fingers growing more quickly along the boundary as predicted by Yang and Yortsos.24 

Four fingering regimes have been identified and studied, namely,

  1. Very early time, when longitudinal diffusion dominates over finger growth.

  2. Early time, when fingers start to grow. This is the flow regime typically investigated using perturbation analysis.

  3. Intermediate time, when multiple fingers grow and interact. This regime can only be investigated using numerical simulation or physical experiments because of the nonlinear behavior.

  4. Late time, when all the fingers, but one, have faded or coalesced.

The number of fingers predicted by the simulations at early time was approximately half the number expected from the analytical solutions obtained from perturbation analysis assuming a sharp interface initially. In addition, the growth rate of those fingers was much lower than predicted by this linear stability analysis. This is consistent with the numerical solutions to the perturbation analysis presented by Tan and Homsy8 for the case when the initially sharp interface spreads due to longitudinal diffusion. Our analysis of their numerical solutions of this linear stability problem showed that the number of initial fingers and the growth rate of those fingers were reduced significantly in this case. Nonetheless, the linear dependence of the number of fingers on PeT and their variation with diffusion anisotropy agreed with the predictions of the perturbation analyses assuming a sharp initial front.

At intermediate and late times, we have shown that the number of fingers over time scale with the most dangerous wavenumber and the growth rate of that wavenumber determined from early time stability analysis of a sharp front presented by Tan and Homsy.8 This is despite the fact that the analysis assumes linear behavior and does not predict the correct number of initial fingers or their growth rate when the front is smearing with longitudinal diffusion. Initially, the number of fingers reduces as 1/t, but at later times, they decay according to 1/t. The along flow spreading of the fingers scales linearly with time in these regimes. This suggests that the empirical models based on modified fractional flow, such as the Koval25 and Todd and Longstaff26 models, can be used to model the average flow behavior in field scale simulation of enhanced oil recovery of CO2 sequestration. We have also shown that the time at which all the fingers merge or fade into one is independent of the viscosity ratio. It depends only on the aspect ratio of the system and the transverse Peclet number.

Further work is needed to explore the impact of velocity dependent dispersion on these findings as well as to explore the fingering scaling and dynamics in radial flows. It would also be useful to design and perform a series of well characterized experiments in pseudo 2D bead packs where the details of the fingering patterns are recorded from very early time and compare these quantitatively with simulations. The systematic investigation of 3D viscous fingering is also an area that would merit further investigation although we note that using numerical simulation to do this remains challenging due to the fine grid resolution required and thus the central processing unit (CPU) time required to perform the simulations. We are investigating the application of advanced numerical methods such as dynamic adaptive mesh optimization in conjunction with parallel processing to achieve this.54 

1.
P. G.
Saffman
and
G.
Taylor
, “
The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid
,”
Proc. R. Soc. London, Ser. A
245
,
312
329
(
1958
).
2.
G.
Rousseaux
,
A.
De Wit
, and
M.
Martin
, “
Viscous fingering in packed chromatographic columns: Linear stability analysis
,”
J. Chromatogr. A
1149
,
254
273
(
2007
).
3.
A.
Kopp
,
H.
Class
, and
R.
Helmig
, “
Investigations on CO2 storage capacity in saline aquifers. Part 1. Dimensional analysis of flow processes and reservoir characteristics
,”
Int. J. Greenhouse Gas Control
3
,
263
276
(
2009
).
4.
M. G.
Gerritsen
and
L. J.
Durlofsky
, “
Modeling fluid flow in oil reservoirs
,”
Annu. Rev. Fluid Mech.
37
,
211
238
(
2005
).
5.
J. P.
Ennis-King
and
L.
Paterson
, “
Role of convective mixing in the long-term storage of carbon dioxide in deep saline formations
,”
SPE J.
10
,
349
356
(
2006
).
6.
H. E.
Huppert
and
J. A.
Neufeld
, “
The fluid mechanics of carbon dioxide sequestration
,”
Annu. Rev. Fluid Mech.
46
,
255
272
(
2013
).
7.
C. T.
Tan
and
G. M.
Homsy
, “
Simulation of nonlinear viscous fingering in miscible displacement
,”
Phys. Fluids
31
,
1330
1338
(
1988
).
8.
C. T.
Tan
and
G. M.
Homsy
, “
Stability of miscible displacements in porous media: Rectilinear flow
,”
Phys. Fluids
29
,
3549
3556
(
1986
).
9.
A.
De Wit
and
G. M.
Homsy
, “
Viscous fingering in periodically heterogeneous porous media. I. Formulation and linear instability
,”
J. Chem. Phys.
107
,
9609
9618
(
1997
).
10.
A.
De Wit
,
Y.
Bertho
, and
M.
Martin
, “
Viscous fingering of miscible slices
,”
Phys. Fluids
17
,
1
9
(
2005
), e-print arXiv:0508080 [physics].
11.
W. B.
Zimmerman
and
G.
Homsy
, “
Viscous fingering in miscible displacements: Unification of effects of viscosity contrast, anisotropic dispersion, and velocity dependence of dispersion on nonlinear finger propagation
,”
Phys. Fluids A
4
,
2348
2359
(
1992
).
12.
M. A.
Christie
and
D. J.
Bond
, “
Detailed simulation of unstable processes in miscible flooding
,”
SPE Reservoir Eng.
2
,
514
522
(
1987
).
13.
C. Y.
Chen
and
E.
Meiburg
, “
Miscible porous media displacements in the quarter five-spot configuration. Part 1. The homogeneous case
,”
J. Fluid Mech.
371
,
233
268
(
1998
).
14.
S. A.
Abdul Hamid
,
A.
Adam
,
M. D.
Jackson
, and
A. H.
Muggeridge
, “
Impact of truncation error and numerical scheme on the simulation of the early time growth of viscous fingering
,”
Int. J. Numer. Methods Fluids
89
,
1
15
(
2019
).
15.
M. N.
Islam
and
J.
Azaiez
, “
Fully implicit finite difference pseudo-spectral method for simulating high mobility-ratio miscible displacements
,”
Int. J. Numer. Methods Fluids
47
,
161
183
(
2005
).
16.
G.
Scovazzi
,
M. F.
Wheeler
,
A.
Mikelic
, and
S.
Lee
, “
Analytical and variational numerical methods for unstable miscible displacement flows in porous media
,”
J. Comput. Phys.
335
,
444
496
(
2017
).
17.
A.
Adam
,
D.
Pavlidis
,
J. R.
Percival
,
P.
Salinas
,
I. C.
London
,
R. D.
Loubens
,
C. C.
Pain
,
A. H.
Muggeridge
, and
M. D.
Jackson
, “
Dynamic mesh adaptivity for immiscible viscous fingering
,” in
SPE Reservoir Simulation Conference
,
2017
.
18.
R. L.
Chuoke
,
P.
van Meurs
, and
C.
van der Poel
, “
The instability of slow, immiscible, viscous liquid-liquid displacements in permeable media
,”
Pet. Trans., AIME
216
,
188
194
(
1959
).
19.
R. L.
Perrine
, “
Stability theory and its use to optimize solvent recovery of oil
,”
Soc. Pet. Eng. J.
1
,
9
16
(
1961
).
20.
J. W.
Gardner
and
J. G. J.
Ypma
, “
An investigation of phase behavior-macroscopic bypassing interaction in CO2 flooding
,”
Soc. Pet. Eng. J.
24
,
508
520
(
1984
).
21.
W. B.
Zimmerman
and
G. M.
Homsy
, “
Nonlinear viscous fingering in miscible displacement with anisotropic dispersion
,”
Phys. Fluids A
3
,
1859
1872
(
1991
).
22.
J. S.
Nijjer
,
D. R.
Hewitt
, and
J. A.
Neufeld
, “
The dynamics of miscible viscous fingering from onset to shutdown
,”
J. Fluid Mech.
837
,
520
545
(
2018
).
23.
M. I.
Morris
and
R. C.
Ball
, “
Renormalization of miscible flow functions
,”
J. Phys. A
23
,
4199
4209
(
1990
).
24.
Z.
Yang
and
Y. C.
Yortsos
, “
Effect of no-flow boundaries on viscous fingering in porous media of large aspect ratio
,”
Soc. Pet. Eng. J.
3
,
285
292
(
1998
).
25.
E.
Koval
, “
A method for predicting the performance of unstable miscible displacement in heterogeneous media
,”
Soc. Pet. Eng. J.
3
,
145
154
(
1963
).
26.
M.
Todd
and
W.
Longstaff
, “
The development, testing and application of a numerical simulator for predicting miscible flood performance
,”
J. Pet. Technol.
24
,
874
882
(
1972
).
27.
F.
Fayers
,
M.
Blunt
, and
M.
Christie
, “
Comparison of empirical viscous-fingering models and their calibration for heterogeneous problems
,”
SPE Reservoir Eng.
7
,
195
203
(
1992
).
28.
K. S.
Sorbie
,
H. R.
Zhang
, and
N. B.
Tsibuklis
, “
Linear viscous fingering: New experimental results, direct simulation and the evaluation of averaged models
,”
Chem. Eng. Sci.
50
,
601
616
(
1995
).
29.
Z. M.
Yang
,
Y. C.
Yortsos
, and
D.
Salin
, “
Asymptotic regimes in unstable miscible displacements in random porous media
,”
Adv. Water Resour.
25
,
885
898
(
2002
).
30.
S.
Malhotra
,
M. M.
Sharma
, and
E. R.
Lehman
, “
Experimental study of the growth of mixing zone in miscible viscous fingering
,”
Phys. Fluids
27
,
014105
(
2015
).
31.
P. M. J.
Tardy
and
J. R. A.
Pearson
, “
A 1D-averaged model for stable and unstable miscible flows in porous media with varying Péclet numbers and aspect ratios
,”
Transp. Porous Media
62
,
205
232
(
2006
).
32.
S. A.
Abdul Hamid
and
A. H.
Muggeridge
, “
Viscous fingering in reservoirs with long aspect ratios
,” in
SPE Improved Oil Recovery Conference
,
Tulsa, Oklahoma, USA
,
2018
.
33.
M. A.
Christie
,
A. H.
Muggeridge
, and
J. J.
Barley
, “
3D simulation of viscous fingering and WAG schemes
,”
SPE Reservoir Eng.
8
,
19
26
(
1993
).
34.
J.
Lohrenz
,
B. G.
Bray
, and
C. R.
Clark
, “
Calculating viscosities of reservoir fluids from their compositions
,”
J. Pet. Technol.
16
,
1171
1176
(
1964
).
35.
E. J.
Peters
,
W. H.
Broman
, and
J. A.
Broman
, “
A stability theory for miscible displacement
,” in
SPE Annual Technical Conference and Exhibition
,
Houston, Texas
,
1984
.
36.
O.
Manickam
and
G. M.
Homsy
, “
Stability of miscible displacements in porous media with nonmonotonic viscosity profiles
,”
Phys. Fluids A
5
,
1356
1367
(
1993
).
37.
R.
Slobod
and
R.
Thomas
, “
Effect of transverse diffusion on fingering in miscible-phase displacement
,”
Soc. Pet. Eng. J.
3
,
9
13
(
1963
).
38.
T. K.
Perkins
,
O. C.
Johnston
, and
N. F.
Hoffman
, “
Mechanics of viscous fingering in miscible systems
,”
Soc. Pet. Eng. J.
5
,
301
317
(
1965
).
39.
J. C.
Bacri
,
N.
Rakotomalala
,
D.
Salin
, and
R.
Wouméni
, “
Miscible viscous fingering: Experiments versus continuum approach
,”
Phys. Fluids A
4
,
1611
1619
(
1992
).
40.
D.
Loggia
,
N.
Rakotomalala
,
D.
Salin
, and
Y. C.
Yortsos
, “
The effect of mobility gradients on viscous instabilities in miscible flows in porous media
,”
Phys. Fluids
11
,
740
742
(
1999
).
41.
H.
Outmans
, “
Nonlinear theory for frontal stability and viscous fingering in porous media
,”
Soc. Pet. Eng. J.
2
,
165
176
(
1962
).
42.
T.
Suekane
,
J.
Ono
,
A.
Hyodo
, and
Y.
Nagatsu
, “
Three-dimensional viscous fingering of miscible fluids in porous media
,”
Phys. Rev. Fluids
2
,
103902
(
2017
).
43.
R.
Blackwell
,
J.
Ryne
, and
W.
Terry
, “
Factors influencing the efficiency of miscible displacement
,”
Pet. Trans., AIME
216
,
1
8
(
1959
).
44.
M. A.
Christie
and
D. J.
Bond
, “
Multidimensional flux-corrected transport for reservoir simulation
,” in
SPE Reservoir Simulation Symposium
,
Bahrain
,
1985
.
45.
S. T.
Zalesak
, “
Fully multidimensional flux-corrected transport algorithms for fluids
,”
J. Comput. Phys.
31
,
335
362
(
1979
).
46.
M. A.
Christie
,
A. D. W.
Jones
, and
A. H.
Muggeridge
, “
Comparison between laboratory experiments and detailed simulations of unstable miscible displacement influenced by gravity
,” in
North Sea Oil and Gas Reservoirs II
, edited by
A. T.
Buller
(
Graham and Trotman
,
London
,
1990
), pp.
245
250
.
47.
G. W.
Davies
,
A. H.
Muggeridge
, and
A. D. W.
Jones
, “
Miscible displacements in a heterogeneous rock: Detailed measurements and accurate predictive simulation
,” in
SPE Annual Technical Conference and Exhibition
,
Dallas
,
1991
.
48.
R. B.
Lantz
, “
Quantitative evaluation of numerical diffusion (truncation error)
,”
Soc. Pet. Eng. J.
11
,
315
320
(
1971
).
49.
T.
Kempers
, “
Dispersive mixing in unstable displacement
,” in
2nd European Conference on the Mathematics of Oil Recovery
, edited by
D.
Guerillot
and
O.
Gospel
(
Editions Technip
,
Paris
,
1990
), Vol. 1, pp.
197
204
.
50.
T. K.
Perkins
and
O. C.
Johnston
, “
A review of diffusion and dispersion in porous media
,”
Soc. Pet. Eng. J.
3
,
70
84
(
1963
).
51.
H. A.
Tchelepi
,
F. M.
Orr
,
N.
Rakotomalala
,
D.
Salin
, and
R.
Woumeni
, “
Dispersion, permeability heterogeneity, and viscous fingering: Acoustic experimental observations and particle-tracking simulations
,”
Phys. Fluids A
5
,
1558
(
1993
).
52.
J. P.
Heller
, “
Onset of instability patterns between miscible fluids in porous media
,”
J. Appl. Phys.
37
(
4
),
1566
1578
(
1966
).
53.
R. A.
Wooding
, “
Growth of fingers at an unstable diffusing interface in a porous medium or Hele-Shaw cell
,”
J. Fluid Mech.
39
,
477
495
(
1969
).
54.
A. E.
Kampsitis
,
A.
Adam
,
P.
Salinas
,
C. C.
Pain
,
A. H.
Muggeridge
, and
M. D.
Jackson
, “
Dynamic adaptive mesh optimisation for immiscible viscous fingering
,”
Comput. Geosci.
(unpublished).