Inertial E ff ects on Dynamics of Immiscible Viscous Fingering in Homogenous Porous Media

We present a comparative study of the onset and propagation dynamics of the fingering phenomenon in uniform porous media with a radial configuration. With the help of the Finite Element Method (FEM)-based 2D simulations and image processing techniques, we investigate finger morphology, growth rate, interfacial length, finger length and the number of fingers which are affected due to inertial forces and convective acceleration in a two-phase porous media flow. We considered a modified Darcy’s law with inertial force coupled with convective acceleration and investigate their impact on interfacial instability with different velocity-viscosity combinations. Interestingly, the consequences of inertial corrections become significant with changes in viscosity at high Reynolds numbers. Due to the intrinsic bifurcation nature of inertial forces in the radial flow geometry, finger morphology is changed mostly at high viscosity ratios. We find that the effects of inertia and convective acceleration are markedly significant at relatively high Reynolds numbers while the interfacial length and the number of fingers—which are important parameters for Enhanced Oil Recovery (EOR)—are most affected by the neglecting of these forces. Moreover, at high Reynolds numbers, the rate of growth of fingering instabilities and the fractal number tend to deviate from that for Darcy’s law.


Introduction
Viscous Fingering (VF) is a form of hydrodynamic instability which evolves when two fluids of unequal viscosities meet such that the less viscous fluid tries to push or penetrate into the more viscous fluid.This phenomenon is also known as the Saffmann-Taylor instability [1].In recent years, the study of immiscible fingering has been investigated due to its wide range of applications in CO 2 sequestration [2,3], Enhanced Oil Recovery (EOR) [4,5], and separation in chromatography [6,7].Depending on the application, sometimes the fingering phenomenon is desirable to enhance process productivity, while sometimes it has a detrimental effect on the final yield.Due to its importance in a wide range of applications related to viscosity-mismatching fluid interactions, many experimental and numerical studies have been conducted to investigate the dynamics and suppression of the fingering phenomenon.Although researchers have focused on various aspects of fingering dynamics and its suppression, it is noticeable that during numerical modelling and their validation with experiments, they have relied on the empirically formulated Darcy's law which was introduced in 1856 [8] and does not take into account inertial effects.While Darcy's equation successfully describes the flow in many applications of porous media, the linear relation between the pressure gradient and the velocity is only valid in specific ranges of fluid velocity.The non-linear relationship between the fluid velocity and the pressure gradient for high speed flows in porous media was first highlighted by Forchheimer in 1901 [9].In addition to the viscous force which was already considered in Darcy's equation, he introduced the inertial force for post-Darcy's zone with a newly-introduced empirical parameter β, where β is a non-Darcy proportionality coefficient which arises in the Forchheimer equation to accommodate inertial effects and has dimension (m −1 ).Before applying the Forchheimer or Darcy's equation, it is important to distinguish between flow regimes for correct applicability of the respective equation.At this point, it is worth mentioning that unlike the pure fluid flow region, researchers have found it difficult to distinguish between the Darcy and post-Darcy zone on the basis of a non-dimensional number.Currently, there exist two different criteria to distinguish between the Darcy and post-Darcy zone.The first criterion-which is similar to non-porous media flow-is related to the evaluation of the Reynolds number and is considered to identify post-Darcy regime.To evaluate the Reynolds number, different studies define it in different ways and set a criterion for the post-Darcy region accordingly [10,11].More recently, the second criterion for recognizing the post-Darcy zone was presented by Zeng et al. [12].The criterion to consider inertial effects in porous media flow is given by dimensionless Forchheimer number, such that Fo = kβν µ ≥ 0.11.Although the value of the non-Darcy coefficient β is obtained through experiments, its nature in terms of variability for different systems has been widely discussed in various studies.Earlier it was considered that for a specific system, β remains constant for all flow rates.However, Baree et al. [13] showed that for the Forchheimer equation, β changes with flow speeds and does not remain constant, thus making it more challenging to incorporate inertial effects in porous media flow.Apart from the Forchheimer equation, several other models have been presented to incorporate the effect of inertial forces in porous media flow.One such model is the modified Hagen-Poiseuille model [14].It is a macroscale model and is based on the hydraulic radius definition.However, just like previous models, it is also based on empirical constants, and the ultimate expression is entirely based on empirical results.After modifications by Burke Plummer [15], and then Ergun et al. [10], the flow equation took the form as: where A and B are empirical constants.It can be observed that the inertial term appearing in Ergun's equation is non-linear, however-despite the non-linearity-it remains an empirical model with limited application for the packed-bed columns.So far, in order to incorporate inertial effects, all the above models can be classified as empirical models.Although these models have integrated inertial terms based on experiments, they lack the theoretical basis in their genesis.At this point, it should also be emphasized that in porous media flow, the instabilities that arise due to inertia may be lesser in magnitude due to small flow rates, but can travel very quickly and grow up in the flow field, and cause the fingering phenomenon to take an effective role [16].To take into account the non-linear convective term along with inertial effects in porous media flow, we found two models which were developed purely on theoretical basis.The first theoretical model was presented by Rabaud et al. [17].They performed different experiments for immiscible displacement and found that the interface was destabilized after a critical value of gas flow.They also observed the propagation and consequent growth of waves in the flow direction with instabilities.After a linear stability analysis, they added a viscous term to the Euler equation and introduced the new equation as: The second theoretical model was proposed by Quil [18].It was purely a theoretical model where he started with the complete Navier-Stokes equation in three dimensions and introduced inertial It should be noted that in the above two equations, different coefficients appear before the non-linear convective terms.Although, both models are derived from averaging the Navier-Stokes equation across the gap of the Hele-Shaw cell, however, both models differ in their averaging techniques.Rabaud et al. used the method of weighted residuals for gap averaging with single test function, while Quil reached to the different coefficients with multiple test functions and better convergence.Nevertheless, it can be seen that the above two equations take into account the inertial as well as the convective terms to model flow in porous media and are free of any empirical constants.Although these equations were derived for a single-phase flow, they may be extended to multiphase flow by using the same analogy of relative permeability as is used for Darcy's law.
Several studies were performed to study the dynamics of the fingering phenomena with and without inertial effects.More recently to incorporate inertial effects, Miranda et al. developed computational models to study the effects of inertia on fingering [19,20].For a rectangular Hele-Shaw cell, they used Darcy's law modified by Gonret [19] and found that inertia had a stabilizing role for the linear regime and it tends to widen the fingers at the weakly non-linear regime.Even though the vast majority of such studies neglected inertial effects on fingering dynamics, those studies provide valuable insight into the subtleties of viscous fingering.Homsy performed experiments in the Hele-Shaw cell to study the viscous fingering both in miscible and immiscible fluid displacement [21].He comprehensively studied the effect of gravity, density, and viscosity on the onset and propagation of fingering.Extreme dendritic fingering growth was observed when a less dense and less viscous fluid tends to displace a denser and more viscous fluid in a Hele-Shaw cell that was tipped vertically.
Similar to the above studies, several investigations were conducted by researchers to suppress fingering.Though the majority of such studies did not consider inertial effects, they provide crucial information on the finger crescendo.Nase et al. examined the pattern formation of Newtonian liquids for lifting a circular Hele-Shaw cell [22].First, the cell was filled with oil, and then the upper plate was lifted so that air could penetrate from sideways and cause fingering instability such that Saffman-Taylor instability evolution could be captured.This setup allowed them to study the instability in a conserved volume of liquid.They found that the surface forces between the fluids were mainly responsible for the number of fingers propagated in the cell.
Juel et al. [23] replaced one of the plates of the Hele-Shaw cell with an elastic membrane and found that the elastic membrane was able to suppress instabilities.The resulting Fluid Structure Interaction (FSI) was found to be responsible for changing the interfacial patterns, which caused a delay in the onset of fingering.Further on FSI, Tahmasebi et al. conducted comprehensive computational studies for multiphase flow in porous media and found that multiphase flow characteristics are modified when stress is applied to the solid porous matrix which can alter dynamics of fluid-fluid interaction [24][25][26].Through various Finite Element Method (FEM)-based models, they found that morphology of solid matrix becomes very crucial during confining stress condition.
Another effort in suppressing instabilities using flow geometry was made by Stone et al. [27].They investigated the effect of tapered geometry on the onset and propagation of fingering and found that a small taper was helpful in suppressing the fingering phenomenon owing to the new parameter in the dispersion relation.The parameter identified was the ratio between the depth gradient to the capillary number.This geometrically mediated control of instability was demonstrated for a classical model system where gas was injected at a constant rate and displaced the viscous wetting fluid.In another study, Stone et al. employed a time-dependent permeability field [28] for finger suppression.During the fluid flow, they lifted the upper plate of a radial Hele-Shaw cell and found that it was also effective in suppressing the fingering phenomenon.They identified that by choosing the right strategy to change the gap thickness in the radial Hele-Shaw cell, fingering could either be suppressed entirely or tip-splitting fingers could be reduced to non-splitting fingers.
Another interesting study to suppress fingering was reported by Ruben et al. [29,30] where they found a change in the fluid-fluid interaction phenomenon by changing the wetting properties of porous media.Through the chemical treatment of the solid surface, they systematically altered the wettability properties in terms of contact angle and found that higher contact angles tend to bring stability to fluid-displacement.
Although numerical models were presented to suppress the fingering phenomenon, all these numerical strategies lack the effects of inertial forces on finger suppression.Although for many porous media flow applications, the assumption of negligible inertial effects remains valid but owing to the relatively high-speed flows in some subsurface flow applications, inertial effects become significant.These applications include costly operations of oil/gas flow near well-bore, liquid waste injection [31], and flow through packed bed reactors which are used in downstream operations [32].In this work, we include inertial corrections along with the convective acceleration in porous media flow and present comparative analysis of viscous fingering phenomenon for immiscible flows.Instead of conventional empirical models, we use a theoretical model to investigate hydrodynamic instabilities arising due to inertial forces.We formulate a 2D numerical model, both with and without inertial effects and investigate the effect of different Reynolds numbers on different aspects of the viscous fingering phenomenon.

Mathematical Modelling
To study the effect of inertial forces on viscous fingering, we consider a 2D Hele-Shaw cell geometry for radial flow such that the outer radius of the geometry is 25 times the inner inlet radius (see Figure 1.) cell, fingering could either be suppressed entirely or tip-splitting fingers could be reduced to non-splitting fingers.
Another interesting study to suppress fingering was reported by Ruben et al. [29,30] where they found a change in the fluid-fluid interaction phenomenon by changing the wetting properties of porous media.Through the chemical treatment of the solid surface, they systematically altered the wettability properties in terms of contact angle and found that higher contact angles tend to bring stability to fluid-displacement.
Although numerical models were presented to suppress the fingering phenomenon, all these numerical strategies lack the effects of inertial forces on finger suppression.Although for many porous media flow applications, the assumption of negligible inertial effects remains valid but owing to the relatively high-speed flows in some subsurface flow applications, inertial effects become significant.These applications include costly operations of oil/gas flow near well-bore, liquid waste injection [31], and flow through packed bed reactors which are used in downstream operations [32].In this work, we include inertial corrections along with the convective acceleration in porous media flow and present comparative analysis of viscous fingering phenomenon for immiscible flows.Instead of conventional empirical models, we use a theoretical model to investigate hydrodynamic instabilities arising due to inertial forces.We formulate a 2D numerical model, both with and without inertial effects and investigate the effect of different Reynolds numbers on different aspects of the viscous fingering phenomenon.

Mathematical Modelling
To study the effect of inertial forces on viscous fingering, we consider a 2D Hele-Shaw cell geometry for radial flow such that the outer radius of the geometry is 25 times the inner inlet radius (see Figure 1.) Fluid I with dynamic viscosity µ1 and density ρ1 enters through the inlet and displaces Fluid II which has dynamic viscosity and density of µ2 and ρ2, respectively.It is assumed that µ1 < µ2.Considering 2D flow in a Hele-Shaw cell, the system of equations for 2-phase flow can be given by, where, ϕ, S and u indicates porosity, fluid saturation, and velocity, respectively.Saturation of both fluids is conserved by the equation: Fluid I with dynamic viscosity µ 1 and density ρ 1 enters through the inlet and displaces Fluid II which has dynamic viscosity and density of µ 2 and ρ 2 , respectively.It is assumed that µ 1 < µ 2 .Considering 2D flow in a Hele-Shaw cell, the system of equations for 2-phase flow can be given by, Fluids 2019, 4, 79 where, φ, S and u indicates porosity, fluid saturation, and velocity, respectively.Saturation of both fluids is conserved by the equation: If we don't consider the inertial and convective acceleration effects, Darcy's model for velocity is given as, In the above equation, p is pressure, k is absolute permeability of medium and k ri is relative permeability.For simplicity, we consider relative permeability as k ri = S 2 i in this work.However, when we consider the inertial corrections, instead of using Darcy's model, we use Quil's model (called the modified Darcy's model in current work), which in the case of immiscible multiphase flow is given by, 6 5 In Equation ( 4), along with other terms, convective acceleration term → u i is also included which manifests itself with inertial effects in the modified Darcy's law.Therefore, the modified Darcy's model incorporates both the effects of inertial forces and convective acceleration.So, we formulated two numerical models: One with Darcy velocity without inertial effects (Equation (1-3)), and the other with the modified Darcy's model encompassing the combined effects of inertia and convective acceleration (Equation ( 1), ( 2) and ( 4)).We generate both models in FEM based COMSOL Multiphysics 5.3 using the 'Coefficient Form PDE" module and performed transient numerical computations.It can also be observed that in addition to convective acceleration, Equation (4) also contains additional terms of time derivative which-combined with convective acceleration-needs discretization similar to the Navier-Stokes equation.This ultimately demands high computing power in terms of speed and memory.As reported by Raju et al. [33], we used Multifrontal Massively Parallel sparse direct Solver (MUMPS) and implemented the self-adaptive time stepping method to solve these equations.
Moreover, to solve the above systems of partial differential equations, it is crucial to choose appropriate initial and boundary conditions.Initially, at time t = 0, we consider only Fluid II inside the flow domain in the stationary state.Fluid I enters the flow domain from inlet such that boundary conditions for velocity and saturation at r = r i are: (i) S 2 = 1 and (ii) → n is the unit vector perpendicular to inlet surface.Furthermore, we impose the boundary conditions of steady pressure and zero saturation flux at the outlet such that the boundary conditions at r = r o are given as: (iii) p = 0 and (iv) → n .∇Si = 0. Before analyzing the results, we find it apposite to discuss a comparison of our numerical model with the experiment and its mesh sensitivity.We set parameters used in the experiment by Ramachandran [34] in our numerical model and compared the results with experimental observations (see Figure 2).We started with a coarse mesh resolution and found that no instabilities were appearing with poor grid resolution.However, as we increased the mesh resolution, the fingering phenomenon starts appearing until we found a close comparison between the experimental result and numerical simulation with an extremely fine mesh resolution (see Figure 3).This shows that the fingering phenomenon is sensitive to the mesh resolution.However, it should be noted here that this is not a new finding and separate independent mesh sensitivity analysis have already been reported in literature [4,35,36], and our finding confirms the mesh dependence phenomenon.We base all our results in this work on extremely fine mesh resolution.Ramachandran [34].

Results and Discussion
In this section, we present our results for the noticeable features that we observed during the comparative study of the dynamics of fingers morphology inferred by the modified Darcy's law.It is already known that the flow in porous media has very low velocities such that their Reynolds number is of the order, 10 −4 -10 −6 [37].However, inertial effects are more prominent at higher Reynolds numbers [19]; therefore, for the current work we perform numerical simulations for Reynolds number of orders of 10  Ramachandran [34].

Results and Discussion
In this section, we present our results for the noticeable features that we observed during the comparative study of the dynamics of fingers morphology inferred by the modified Darcy's law.It is already known that the flow in porous media has very low velocities such that their Reynolds number is of the order, 10 −4 -10 −6 [37].However, inertial effects are more prominent at higher Reynolds numbers [19]; therefore, for the current work we perform numerical simulations for Reynolds number of orders of 10

Results and Discussion
In this section, we present our results for the noticeable features that we observed during the comparative study of the dynamics of fingers morphology inferred by the modified Darcy's law.It is already known that the flow in porous media has very low velocities such that their Reynolds number is of the order, 10 −4 -10 −6 [37].However, inertial effects are more prominent at higher Reynolds numbers [19]; therefore, for the current work we perform numerical simulations for Reynolds number of orders of 10 the plates and denotes characteristic length.Moreover, in order to generalize our results, we define another parameter viscosity ratio, R defined as R = µ 2 µ 1 , and present our results with a dimensionless time parameter T = U inlet t h . In this study, the fingering phenomenon is studied for values of R ranging from 10-1000.We present our results in two sections: The first section relates to qualitative analysis where we compare the evolution of fingers with and without inertial and convection acceleration effects alongside various combinations of Re and R. In the second section, we quantify effects of the modified Darcy's law and present comparative analysis for Darcy's and the modified Darcy's models.

Qualitative Analysis
First, we compare saturation contour plots of Fluid I for Darcy's and the modified Darcy's models.After performing simulations with various velocity-viscosity combinations, we find that the saturation of fingering travels inside the flow domain in various levels.Levels with low saturation of fluid travel with high speeds while high saturation levels travel at relatively low speed and remain far from the displaced fluid (see Figure 4).Moreover, fingers morphology also changes with the saturation level.For example, it can be seen that the high saturation level fluid stays near the flow inlet with uniform spread, while the low saturation level fluid travels outwards quickly with penetrating effect on Fluid II.This fact remains true both for Darcy's and the modified Darcy models.However, at the moderate Reynolds number and high viscosity ratio (Re = 0.0219 and R = 1000 in this case), more finger split is evident for the modified Darcy's model.Though for rectilinear geometric flow configurations [19], inertia has a stabilizing effect, but bifurcation of fingering is observed.
'h' is gap thickness between the plates and denotes characteristic length.Moreover, in order to generalize our results, we define another parameter viscosity ratio, R defined as ( ) , and present our results with a dimensionless time parameter ( ) In this study, the fingering phenomenon is studied for values of R ranging from 10-1000.We present our results in two sections: The first section relates to qualitative analysis where we compare the evolution of fingers with and without inertial and convection acceleration effects alongside various combinations of Re and R. In the second section, we quantify effects of the modified Darcy's law and present comparative analysis for Darcy's and the modified Darcy's models.

Qualitative Analysis
First, we compare saturation contour plots of Fluid I for Darcy's and the modified Darcy's models.After performing simulations with various velocity-viscosity combinations, we find that the saturation of fingering travels inside the flow domain in various levels.Levels with low saturation of fluid travel with high speeds while high saturation levels travel at relatively low speed and remain far from the displaced fluid (see Figure 4).Moreover, fingers morphology also changes with the saturation level.For example, it can be seen that the high saturation level fluid stays near the flow inlet with uniform spread, while the low saturation level fluid travels outwards quickly with penetrating effect on Fluid II.This fact remains true both for Darcy's and the modified Darcy models.However, at the moderate Reynolds number and high viscosity ratio ( 0.0219 Re = and 1000 R = in this case), more finger split is evident for the modified Darcy's model.Though for rectilinear geometric flow configurations [19], inertia has a stabilizing effect, but bifurcation of fingering is observed.It is already known from fluid mobility analysis that less viscous fluids travel much faster than more viscous fluids.In our case, we also observe the same phenomenon (see Figure 5), but we find that this phenomenon takes place at the expense of stability.More viscous fluids travel with a highly stabilized fluid-fluid interface such that the finger-spread area remains more uniform and less finger-split occurs.On the other hand, less viscous fluids are not stable and tend to produce a larger size of fingering pattern with more finger-splits.This observation remains unchanged when the inertial effect and convective acceleration are considered in the modified Darcy's model.However, the phenomenon of increased stability with a viscous fluid is strengthened due to Darcy's law modification.It is evident from Figure 5 that at 100 R = , inertia and convective acceleration widen the finger size whereas at 1000 R = negligible difference can be noticed as an impact of Darcy's law modification.
It is already known from fluid mobility analysis that less viscous fluids travel much faster than more viscous fluids.In our case, we also observe the same phenomenon (see Figure 5), but we find that this phenomenon takes place at the expense of stability.More viscous fluids travel with a highly stabilized fluid-fluid interface such that the finger-spread area remains more uniform and less finger-split occurs.On the other hand, less viscous fluids are not stable and tend to produce a larger size of fingering pattern with more finger-splits.This observation remains unchanged when the inertial effect and convective acceleration are considered in the modified Darcy's model.However, the phenomenon of increased stability with a viscous fluid is strengthened due to Darcy's law modification.It is evident from Figure 5 that at R = 100, inertia and convective acceleration widen the finger size whereas at R = 1000 negligible difference can be noticed as an impact of Darcy's law modification.Evolution of fingers changes with viscosity.Figure 6 shows the evolution of the fluid interface with different time instants.We note that at some locations, fingering with the modified Darcy's model tend to split more than fingering with the Darcy's model, but this trend is inconsistent and hence cannot be generalized.However, as we decrease viscosity while keeping inlet velocity constant, we observe that fingers are more unstable, both with and without Darcy's law modification (see Figure 7).Evolution of fingers changes with viscosity.Figure 6 shows the evolution of the fluid interface with different time instants.We note that at some locations, fingering with the modified Darcy's model tend to split more than fingering with the Darcy's model, but this trend is inconsistent and hence cannot be generalized.However, as we decrease viscosity while keeping inlet velocity constant, we observe that fingers are more unstable, both with and without Darcy's law modification (see Figure 7).Evolution of fingers changes with viscosity.Figure 6 shows the evolution of the fluid interface with different time instants.We note that at some locations, fingering with the modified Darcy's model tend to split more than fingering with the Darcy's model, but this trend is inconsistent and hence cannot be generalized.However, as we decrease viscosity while keeping inlet velocity constant, we observe that fingers are more unstable, both with and without Darcy's law modification (see Figure 7).We also compared the evolution of both models in different frames with low Fluid II viscosity and relatively high velocity (see Figure 8.)Although there were fewer numbers of finger splits for  We also compared the evolution of both models in different frames with low Fluid II viscosity and relatively high velocity (see Figure 8.)Although there were fewer numbers of finger splits for the modified Darcy's model, no significant changes were observed.We also compared the evolution of both models in different frames with low Fluid II viscosity and relatively high velocity (see Figure 8.)Although there were fewer numbers of finger splits for the modified Darcy's model, no significant changes were observed.

Quantitative Analysis
In this section, we enumerate the inertial effects of inertia and convective acceleration on fingering in terms of fractal number, interfacial length, number of fingers, growth rate, and fingers length.For ease of view and comparison purpose, we present our results in order in Figure 9 for the Reynolds number to make comparative analysis easier.
The fractal number is a measure of a shape's ability to repeat itself with respect to its change in scale.Non-dimensional fractal number (Df) provides a meticulous capability to study morphology and pattern formations of the fingering phenomenon [38,39].In the current work, quantification of fingers morphology and pattern formation was carried out using dimensionless fractal number analysis.As part of the post-processing route, we acquired color images of fingering patterns from COMSOL and processed them to convert to 8-bit binary images.The fractal dimensions of binary images were computed by the box-counting method using FracLac (a plugin of a freeware program ImageJ [40]).Overall, we observed that an increase in the Reynolds number increases the fractal number, irrespective of inclusion of inertial and convective acceleration effects.Moreover, the fractal

Quantitative Analysis
In this section, we enumerate the inertial effects of inertia and convective acceleration on fingering in terms of fractal number, interfacial length, number of fingers, growth rate, and fingers length.For ease of view and comparison purpose, we present our results in order in Figure 9 for the Reynolds number to make comparative analysis easier.
The fractal number is a measure of a shape's ability to repeat itself with respect to its change in scale.Non-dimensional fractal number (D f ) provides a meticulous capability to study morphology and pattern formations of the fingering phenomenon [38,39].In the current work, quantification of fingers morphology and pattern formation was carried out using dimensionless fractal number analysis.As part of the post-processing route, we acquired color images of fingering patterns from COMSOL and processed them to convert to 8-bit binary images.The fractal dimensions of binary images were computed by the box-counting method using FracLac (a plugin of a freeware program ImageJ [40]).Overall, we observed that an increase in the Reynolds number increases the fractal number, irrespective of inclusion of inertial and convective acceleration effects.Moreover, the fractal number also increases with the evolution of fingers.This conforms the results by Pons et al. [38].However, it is evident that the modification in Darcy's law tends to decrease the fractal number at the high Reynolds number.Since high values of (D f ) signify the capability of fingers to split into similar shapes, we observe that inertial forces combined with convective acceleration tend to produce fingers of shape which are less similar to the parent fingers.Similar findings without quantification were reported by Miranda et al.where they refer to these non-similar fingers as 'interfacial lobes' branching out sideward [41].
Another critical parameter which plays a vital role in reactive flows is the interfacial length.Again, we used image processing freeware ImageJ to compute the length interface in terms of the number of pixels and compared them to different cases (see Figure 9b).It is interesting to note that at low Reynolds numbers, non-consideration of inertial effects and convective acceleration in Darcy's model underestimates interfacial length.Clearly, it can be observed that with the propagation of fingers, the ratio of interfacial length with and without Darcy's law-modification is increased.This can have considerable effects on the estimation of Enhanced Oil Recovery using reactive polymer flows.Such flows are governed by low Reynolds numbers of order 10 -5 where contact area plays a vital role in total oil recovery.In such cases, application of Darcy's law can potentially misjudge the total oil recovery.However, it is evident that the modification in Darcy's law tends to decrease the fractal number at the high Reynolds number.Since high values of (Df) signify the capability of fingers to split into similar shapes, we observe that inertial forces combined with convective acceleration tend to produce fingers of shape which are less similar to the parent fingers.Similar findings without quantification were reported by Miranda et al.where they refer to these non-similar fingers as 'interfacial lobes' branching out sideward [41].For most applications, it is desirable to have reduced growth rate, fingers length and number of fingers.To get these parameters from our numerical simulations, we first converted interfacial contours to Cartesian coordinates and then processed the plot by employing several algorithms to determine the mean length, quantity and growth rates of fingers.It is evident that at low Reynolds numbers, the number of bifurcations increases with consideration of inertial effects and convective acceleration (see Figure 9c).This also explains the increase in interfacial length at low Reynolds numbers.However, as we increase Reynolds numbers, the number of fingers is increased to approximately twofold, with and without Darcy's law modification.From Figure 9c1,c2, it is difficult to extract an accurate trend between the Darcy's and the modified Darcy's law using one-on-one comparison.This arises due to the increased noise in signals with increased Reynolds numbers where we found it difficult to differentiate between minor fingers and noise.However, by using the trend of interfacial length, we can deduce that the number of fingers produced is actually no lesser at high Reynolds numbers, but obviously, fingers are not produced at the same rate as in the case of low Reynolds numbers.
Inertial forces coupled with convective acceleration also tend to decrease the growth rate of fingers, but these effects become more prominent in later stages of fingers evolution (see Figure 9d,e).As we can see that growth rate is always negative, which implies that wavenumber (n) of the system has already passed the most dangerous wavenumber (n max ).However, we can see that the effect of growth rate becomes significant only at high Reynolds numbers.Moreover, the length of fingers also tends to be affected at later stages of fingers propagation.Though we observe shorter fingers with inertial effects and convective acceleration for all Reynolds numbers considered in this study, yet we can see that the extent of fingers length grows shorter with an increase in Reynolds numbers.This observation has its implications in EOR using water flooding where the length of fingers affects oil recovery.However, for the fact that EOR takes place at low Reynolds numbers during water flooding, fingers lengths are not expected to change much in this case.

Conclusions
In this study, we attempted to find the effects of inertia and convective acceleration on fingers morphology in homogenous porous media for radial configuration.We selected a theoretical model to consider the inertial effects in the modified Darcy's model and compared both Darcy's and the modified Darcy's models for different combinations of viscosity and inlet velocity of the displacement fluid.The salient features of our study are:

•
The trend of saturation distribution remains the same, both for Darcy's and the modified Darcy's law with inertial corrections and convective acceleration.High saturation levels of invading fluid stay near the inlet while levels with low saturation travel far from the inlet.

•
While considering the modification in Darcy's law, the fractal number is found to decrease as the Reynolds number increases.

•
The interfacial area is observed to be underestimated by Darcy's law, especially at low Reynold's numbers.

•
Modification in Darcy's law affects fingers morphology.At low viscosity ratios, fingers tend to widen while this effect is relatively small at high viscosity ratios.

•
Length and growth of fingers are affected mostly at later stages of the fingers evolution, while the growth rate of fingers decreases with an increase in inertial effects and convective acceleration.
Although in our numerical study, we considered several combinations of Reynolds numbers and viscosity ratios, due to associated computational cost, results were not produced for a larger extent of these combinations with time.It may be interesting to perform a non-dimensional analytical study with inertial terms and investigate the effect of inertial terms with different combinations.The results of the analytical study-coupled with the extraction of important parameters through image processing techniques-should provide more general consequences of the modified Darcy's law.

Fluids
perturbative fashion.The modified Darcy's equation for porous media flow was given as: 6 5

Figure 1 .
Figure 1.Schematic of radial flow domain where less viscous fluid (red) displaces more viscous fluid (blue).

Figure 1 .
Figure 1.Schematic of radial flow domain where less viscous fluid (red) displaces more viscous fluid (blue).

Figure 3 .
Figure 3. Mesh sensitivity analysis for comparison of numerical model with the experiment by

Figure 3 .
Figure 3. Mesh sensitivity analysis for comparison of numerical model with the experiment by

Figure 3 .
Figure 3. Mesh sensitivity analysis for comparison of numerical model with the experiment by Ramachandran [34].

Fluids
−1 -10 −3 .The flow domain is considered to have a uniform porosity of φ = 0.5 We change the properties of Fluid I and Fluid II to attain different Reynolds number for a comparative study between Darcy's law and Quil's modified Darcy's law (to be called modified Darcy's law henceforth).It should be emphasized here that the modified Darcy's law encompasses expression of inertial force combined with convective acceleration, hence the results of the modified Darcy's law should be construed in this context.We define Reynold's number as Re = ρUh µ while 'h' is gap thickness between

Figure 4 .
Figure 4. Fingering phenomenon with Fluid I saturation contour plots for (a) Darcy's and (b) the modified Darcy's models for 0.0219 Re = and 1000 R = at 454.5 T = .

Figure 5 .
Figure 5.Comparison of propagation of fingering for fluids with different viscosities for (a) Darcy and (b) the modified Darcy's model for contours of saturation 1 0.02 S = are compared for fluids of viscosity 100 R = (blue) and 1000 R = (red).Both contours are captured for the same inlet velocity at time 454.4T = .

Figure 6 .
Figure 6.Comparison of the evolution of fingering for (a) Darcy's and (b) the modified Darcy's model for contours of saturation 1 0.02 S = for fluids of viscosity ratio 100 R = and Reynolds

Figure 5 .
Figure 5.Comparison of propagation of fingering for fluids with different viscosities for (a) Darcy and (b) the modified Darcy's model for contours of saturation S 1 = 0.02 are compared for fluids of viscosity R = 100 (blue) and R = 1000 (red).Both contours are captured for the same inlet velocity at time T = 454.4.

Fluids 2019, 4 , 13 Figure 5 .
Figure 5.Comparison of propagation of fingering for fluids with different viscosities for (a) Darcy and (b) the modified Darcy's model for contours of saturation 1 0.02 S = are compared for fluids of viscosity 100 R = (blue) and 1000 R = (red).Both contours are captured for the same inlet velocity at time 454.4T = .

Figure 6 .
Figure 6.Comparison of the evolution of fingering for (a) Darcy's and (b) the modified Darcy's model for contours of saturation 1 0.02 S = for fluids of viscosity ratio 100 R = and Reynolds

Fluids 2019, 4 , 13 Figure 7 .
Figure 7.Comparison of the evolution of fingering for (a) Darcy's and (b) the modified Darcy's model contours of saturation 1 0.02 S = for fluids of viscosity ratio 1000 R = and Reynolds

Figure 7 .
Figure 7.Comparison of the evolution of fingering for (a) Darcy's and (b) the modified Darcy's model contours of saturation S 1 = 0.02 for fluids of viscosity ratio R = 1000 and Reynolds number of Re = 2.19 × 10 −3 for time T = 18 (blue), T = 181 (red) and T = 363 (green).

Figure 7 .
Figure 7.Comparison of the evolution of fingering for (a) Darcy's and (b) the modified Darcy's model contours of saturation 1 0.02 S = for fluids of viscosity ratio 1000 R = and Reynolds

Figure 8 .
Figure 8.Comparison of the evolution of Fluid I saturation contour plots for Darcy's (top) and the modified Darcy's model (bottom) for 1000 R = and

Figure 8 .
Figure 8.Comparison of the evolution of Fluid I saturation contour plots for Darcy's (top) and the modified Darcy's model (bottom) for R = 1000 and Re = 2.12 × 10 −2 .

Figure 9 .
Figure 9.Comparison of Darcy's model with the modified Darcy's Model at R = 1000 for (a) fractal number, (b) interfacial length, (c) number of fingers, (d) growth rate and (e) finger length at different Re.