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.

## I. INTRODUCTION

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 (CO_{2}) 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 Ball^{23} and subsequently Yang and Yortsos^{24} 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 Yortsos^{24} 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 Koval^{25} or Todd and Longstaff^{26} 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 *D*_{T}, and longitudinal diffusion *D*_{L} 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.

## II. THEORY

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.

### A. Flow equations

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

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

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

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

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

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

where *D*_{o} 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 Homsy^{21} 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),

### B. Dimensionless scaling

Previous workers have suggested a number of different ways of scaling the system. In particular, we note that Tan and Homsy^{7} and Zimmerman and Homsy^{11} scale all lengths by the longitudinal dispersive length (*D*_{L}/|**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

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

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

The choice of *t*_{D} 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

where *ϵ* = *D*_{T}/*D*_{L}, $PeT=Hvt\varphi DT$ is the transverse Peclet number, $vD=vvt$, and $vt=1H\u222b0Hvxdy$ is the average speed in the direction of flow.

### C. Linear perturbation theory

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

where $\u0101$ is the mean value of the quantity (averaged transverse to the flow) and *a*′ is the fluctuation to the mean such that

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

where $\Phi xD,tD=\Phi xD,0e\omega tD$, i.e., the amplitude of the fluctuation grows at a rate of *ω* and the base state, $\u0101$, is a step function in concentration at *t*_{D} = 0. They also assume periodic boundary conditions.

Tan and Homsy,^{8} amongst others, showed that there is a dimensionless finger wavenumber *k*_{m} with maximum dimensionless growth rate *ω*_{m} and a cutoff wavenumber *k*_{c} 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 Bond^{12}). If the viscosity varies exponentially with concentration [Eq. (4)], then Tan and Homsy^{8} were able to obtain an analytical expression for the initial growth rate (at *t*_{D} = 0) as a function of wavenumber, mobility ratio, and dispersion given by

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

The most dangerous wavenumber *k*_{m} occurs when $d\omega dk=0$. In this case, Eq. (15) can only be solved exactly when diffusion is isotropic, giving

and

Tan and Homsy^{8} provide two approximate solutions for when *ϵ* ≫ 1 (transverse diffusion is much greater than longitudinal),

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

We see that in all cases, the wavenumber scales proportionally to (*RPe*_{T}) and the growth rate of the fingers scales proportionally to (*RPe*_{T})^{2}. Note that in our system of dimensionless parameters, the dimensionless wavenumber is *k* = 2*πn*_{f}, where *n*_{f} is the number of fingers. Manickam and Homsy^{36} showed that the preceding analysis can be extended to any function relating mixture viscosity to concentration by replacing *R* with Λ, where Λ is defined as

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 Thomas^{37}).

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 *k*_{c}, *k*_{m}, 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 Homsy^{8} were able to obtain numerical solutions for the case *R* = 3 and *ϵ* = 1. Figure 3 sketches the overall form of Tan and Homsy’s^{8} results, while Fig. 3(b) shows the number of fingers corresponding to the most dangerous wavenumber (*n*_{f}) as a function of dimensionless advective time, *t*_{D, 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 *t*_{D, A} = 0 and *t*_{D, 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’s^{8} 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 Homsy^{8} are also given. The number of fingers at *t*_{D, A} = 0.04 is approximately half that at *t*_{D, A} = 0, while the growth rate is reduced by a factor of 4.

Pe_{T}, ϵ = 1
. | n_{f}, Eq. (17)
. | n_{f}, t_{D,A} = 0
. | n_{f}, t_{D,A} = 0.04
. | n_{f,0}/n_{f,0.04}
. | ω_{m} Eq. (18)
. | ω_{m} t_{D,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 |

### D. Number of fingers at early time

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 *k*_{m} by

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 *R*^{2} to obtain

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

For any given *ϵ*, there is a linear relationship between $nfR$ and *Pe*_{T}. 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 *ϵ* → ∞,

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.

## III. FINGERING REGIMES

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

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}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.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.

Late time, when a single finger is left and propagates through the system. This was first described by Saffman and Taylor

^{1}and analyzed by Outmans^{41}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 CO_{2} 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.5^{42}) and along the finger for higher viscosity ratios (e.g., *M* = 383^{43}). 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 Homsy^{21} 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 Homsy^{7} 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 Homsy^{11} 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.

## IV. NUMERICAL SIMULATION

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 Bond^{44} 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. Simulation input parameters

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 Homsy^{21}) 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.

. | Very early to . | Intermediate . |
---|---|---|

Parameters . | intermediate times . | to late time . |

Domain aspect ratio, (L/H) | 1 | 10 |

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, Pe_{T} | ||

Ratio of transverse to | 0.2, 1, 5 | 0.25, 0.5, 1, 2 |

longitudinal dispersion, ϵ |

. | Very early to . | Intermediate . |
---|---|---|

Parameters . | intermediate times . | to late time . |

Domain aspect ratio, (L/H) | 1 | 10 |

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, Pe_{T} | ||

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 ∆*x*_{D} = 1/600 = 0.001 67 and ∆*y*_{D} = 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 ∆*x*_{D} = ∆*y*_{D} = 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.

### B. Analysis of results

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, *t*_{D,A} (pore volumes injected for a unit aspect ratio system), given by

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 *L*_{mix},

where *x*_{D}|_{c=0.02} and *x*_{D}|_{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 *c*_{level} 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 *c*_{level} 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. *c*_{level} = 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.

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.

## V. RESULTS

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 *Pe*_{T} = 250, *ϵ* = 1, and *R* = 3. Figure 6(a) shows the number of fingers *n*_{f} as a function of *t*_{D,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.

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 *t*_{D,A} = 0.3, 8 small fingers form. The number of fingers then decays monotonically until at around *t*_{D,A} = 0.9, only 3 fingers remain. Just after *t*_{D,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 *t*_{D,A} = 2.3 when another splitting event occurs, until just before *t*_{D,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 Koval^{25} and Todd and Longstaff^{26} 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.

### A. Very early to early time regime

As noted in Sec. III, at very early time, the sharp front between the displacing and displaced fluids is diffused by longitudinal dispersion *D*_{L}. 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 *Pe*_{T} 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 *Pe*_{T}, *L*_{mix} 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 *R*^{2} [Eq. (30)].

Kempers^{49} argued heuristically that the time for this crossover can be roughly estimated by equating the mixing length induced by longitudinal diffusion,

to the linear spreading based on the Koval^{25} model,

where *M*_{e} is the effective viscosity ratio defined as

Note that we have adapted the Kempers^{49} 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 Koval^{25} 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,

where *l*_{D} 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, *Pe*_{T} = 500, and *ϵ* = 1. Here, *l*_{D} 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 Homsy^{8} at *t*_{D,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 Kempers^{49} 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 *t*_{D,A} ∼ 0.04) than the Kempers^{49} approach, but this is considerably earlier than the transition seen in the simulation (which transitions to early time fingering at around *t*_{D,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 Homsy^{8} for *t*_{D,A} = 0.04 results in a much later estimate of a transition to viscous fingering (at *t*_{D,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 *t*_{D,A}.

Figure 8 shows the number of fingers obtained at early time from simulation, scaled by *R*, as a function of *Pe*_{T} and *ϵ*. These follow the same linear trend as the maximum number of fingers predicted from the perturbation analysis of Tan and Homsy^{8} as shown in Fig. 4. The number of fingers increases with the ratio of transverse to longitudinal dispersion for a given *Pe*_{T} 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.

It can be seen, however, that the gradients of *n*_{f} vs *Pe*_{T} 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 Homsy^{8} 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.

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 Homsy^{21} 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 Heller^{52} 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.

### B. Intermediate time regime

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 *Pe*_{T} = 250, 500, and *ϵ* = 1. The figure shows the concentration distribution at *t*_{D,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 *t*_{D,A}. As expected, initially more fingers start growing sooner and grow more quickly for larger *Pe*_{T} or *R*. The fingering dynamics become nonlinear sooner with more splitting/branching occurring at higher viscosity ratios and *Pe*_{T}.

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

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.

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.

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 Azaiez

^{15}where a pair of fingers adjacent to a dominant, shielding finger bend and merge into its base.Tip splitting (TS), which refers to the instability when a finger-tip splits into two. Tan and Homsy

^{7}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*Pe*_{T}and*ϵ*of the system and the longitudinal concentration gradient at the tip.Side branching (SB) is also observed particularly in the highest mobility ratio, highest

*Pe*_{T}case.

In general, we observe that the frequency of splitting events increases when we have a higher mobility ratio and *Pe*_{T}.

Figure 11 shows the effects of increasing and reducing longitudinal dispersion (*ϵ* = 0.2, 5) when *Pe*_{T} = 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 *Pe*_{T}, 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 *D*_{T}. Although this case is unrealistic (as longitudinal dispersion *D*_{L} is typically higher than *D*_{T} for flow in porous media), it provides further evidence that a smaller *D*_{L} will result in a greater flow instability which is consistent with the early time linear stability analysis.

### C. Scaling of finger number vs time

Figure 12 shows the number of fingers, *n*_{f,} as a function of *t*_{D,A} for *Pe*_{T} = 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 Wooding^{53} showed that the number of fingers reduces diffusively as a function of time, i.e., *n*_{f} ∼ *t*^{−0.5}, and this depended upon *D*_{T}. Although we observed similar behavior for *R* = 3 for 0.1 < *t*_{D,A} < 0.3 [see Fig. 11(b)], for higher mobility ratios, we observed the following:

Initially, the number of fingers decreases approximately as $nf\u223ctD,A\u22121$ [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 $nf\u223cRPeTtD,A\u22121$.After this initial period, for

*R*= 4, 5,*n*_{f}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*Pe*_{T}= 500 case for*t*_{D,A}up to 0.5. This is similar to the observations by Zimmerman and Homsy^{21}who found that they had to increase*Pe*_{T}> 300 before tip splitting occurred in their spectral method simulations.

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

Note that *t*_{D} is the dimensionless diffusive time defined in Eq. (9). Nijjer *et al.*^{22} proposed the same scaling for *n*_{f} but scaled advective time by *Pe*_{T}*R*^{2} 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 *Pe*_{T} = 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.

### D. Late time regime

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 CO_{2} sequestration and enhanced oil recovery.

Figure 14 shows a typical set of results for the case *R* = 3, *Pe*_{T} = 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 *t*_{D,A} = 5 for visual reference. Figures 14(c) and 14(d) show the number of fingers and mixing length as a function of *t*_{D,A}, respectively. Note that the time scale in Fig. 14(d) is logarithmic.

Figure 14(a) shows that the displacement starts with 12 fingers but very quickly (by *t*_{D,A} = 1.8) nine of those fingers have faded due to shielding and coalescence. Two of the remaining fingers experience tip splitting at approximately *t*_{D,A} = 2 and *t*_{D,A} = 2.5, but there are only 2 fingers by *t*_{D,A} = 3. In contradiction to the analysis of Yang and Yortsos^{24} 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 *t*_{D,A} ∼ 0.1. After this time, the number of fingers decays monotonically (except for a brief period when there are tip splitting events between *t*_{D,A} = 2 and *t*_{D,A} = 3) to two fingers at *t*_{D,A} = 3.3. The mixing length grows linearly with *t*_{D,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, *Pe*_{T} = 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 *n*_{f}*/*(*k*_{m}*/*2*π*) against *t*_{D,A}*ω*_{m}, as shown in Fig. 15(a).

For isotropic dispersion, the mixing lengths *L*_{mix} for *Pe*_{T} = 500, 250, 125 exhibit different behaviors in the different fingering regimes, as shown in Fig. 15(b). A lower *Pe*_{T} gives a higher mixing length at very early time (as longitudinal dispersion *D*_{L} is higher) but yields a lower mixing length at intermediate to late times.

Plotting *L*_{mix} against dimensionless advective time for simulations with the same longitudinal dispersion (*D*_{L} = 0.004) but different *Pe*_{T} = 125, 250, 500 in Fig. 15(c) gives identical *L*_{mix} at very early time, as expected, however, at intermediate to late time regimes, *L*_{mix} depends mainly upon *Pe*_{T}, as shown in Fig. 15(d) where we plot the cases with the constant transverse Peclet number (*Pe*_{T} = 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 *Pe*_{T} = 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.

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 suggestion^{7} 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 Homsy^{21} 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 *Pe*_{T}. 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 *t*_{D,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 *D*_{L} results in a steeper concentration gradient in the direction of flow at the front of the finger, but a higher *D*_{T} will tend to spread the finger transverse to flow. There is thus an optimum combination of *D*_{L}, *D*_{T}, and mobility ratio for tip-splitting. In particular, if *D*_{L} 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}

Tan and Homsy^{7} 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\pi km\u221d1tD\omega m$; thus, the time at *k*_{m} which all the fingers have disappear except for 1 can be estimated from

(and there is only one finger left so *n*_{f} = 1). Assuming isotropic diffusion and substituting for *k*_{m} from Eq. (17) and *ω*_{m} from Eq. (18), we obtain

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, *t*_{D,A}, at which the system decays to 1 finger results in

It is interesting that the time does not seem to depend directly upon *R*. Gardner and Ypma^{20} plotted unrecovered oil at *t*_{D,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,

The simulations and experiments were performed in different aspect ratio systems, yet they found that recovery scaled very well with 1/*N*_{T} (see their Fig. 11). They found that the volume of unrecovered oil was constant for 1/*N*_{T} < 0.02 and reduced as *N*_{T} increased for 1/*N*_{T} > 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 *Pe*_{T}, 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 *Pe*_{T}; however, further simulations, supported by physical experiments, would be needed to substantiate this hypothesis.

## VI. CONCLUSIONS

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 Homsy^{8} and Zimmerman and Homsy.^{21} The combination of isotropic diffusion and high Peclet numbers means the finger-tip, for Nijjer *et al.*’s^{22} combination of *R* and *Pe*_{T}, never becomes sufficiently unstable to split. Our study investigated anisotropic diffusion and slightly lower Peclet numbers (*Pe*_{T} ≤ 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,

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

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

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.

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 Homsy^{8} 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 *Pe*_{T} 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 Koval^{25} and Todd and Longstaff^{26} models, can be used to model the average flow behavior in field scale simulation of enhanced oil recovery of CO_{2} 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}

Electronic mail: ann.muggeridge@imperial.ac.uk

## REFERENCES

_{2}storage capacity in saline aquifers. Part 1. Dimensional analysis of flow processes and reservoir characteristics

_{2}flooding