1 Introduction

## Abstract

A hierarchy of mathematical models describing viscosity-stratified flow in a Hele-Shaw cell is constructed. Numerical modelling of jet flow and development of viscous fingers with the influence of inertia and friction is carried out. One-dimensional multi-layer flows are studied. In the framework of three-layer flow the interpretation of the Saffman–Taylor instability is given. A kinematic-wave model of viscous fingering taking into account friction between the fluid layers is proposed. Comparison with calculations on the basis of two-dimensional equations shows that this model allows to determine the velocity of propagation and the thickness of the viscous fingers.

VISCOSITY-STRATIFIED FLOW IN A HELE-SHAW CELL

A. A. Chesnokov, V. Yu. Liapidevskii

[2mm]

Novosibirsk State University,

Pirogova Str. 2, Novosibirsk, 630090, Russia

[2mm] Lavrentyev Institute of Hydrodynamics SB RAS,

Lavrentyev Ave. 15, Novosibirsk, 630090, Russia

e-mails: chesnokov@hydro.nsc.ru, liapid@hydro.nsc.ru

Keywords: Hele-Shaw flow, fingering instability, kinematic-wave model.

## 1 Introduction

A displacement process involving two fluids is often unstable when the displacing fluid has larger viscosity than the displaced one. The resulting instability developing at the interface between two fluids is known as viscous fingering [20, 12]. This instability has received much attention as an archetype of pattern-formation problems and as a limiting factor in the recovery of crude oil. Classical mathematical model describing a Newtonian flow displacement in a Hele-Shaw cell and development of the Saffman–Taylor instability consists of the continuity equation, Darcy’s law and a convection-diffusion equation for the concentration of the displacing fluid [22, 3]. The inertia of fluid may be important for high finger velocities. This leads to the necessity to use more complex nonlinear equations of fluid motion [10]. In the framework of these models instability caused by different velocities of layers movement can be considered. There are several theoretical and experimental studies on the role of inertia in immiscible [7, 8] and miscible [27] displacements. The results reveal that inertia tends to damping viscous fingering. In recent publications [2, 18] the effects of the thickness variation of a Hele-Shaw cell and elasticity of the walls on the process of viscous fingering have been studied. Different types of instability in viscosity-stratified flow have been discussed in [11].

A number of theoretical, numerical, and experimental works devoted to understanding of various aspects of the instabilities is available in the literature [16, 23, 21]. Accurate numerical solution is costly at high Peclet numbers and it is difficult to reproduce the detailed fingering pattern. Numerical simulations [28, 25] show that it is simpler to describe the concentration of solute averaged across the fingers. In such cases, the mixing zone is an important feature to determine the extent of mixing. Despite the considerable work done, the spreading and growth of the mixing zone is an important question that still remains unresolved. Several empirical models are available for the evaluation of mixing zones in unstable, miscible displacements. Two empirical models have been suggested by Koval [13] and Todd and Longstaff [24] to give a basis for computation of miscible displacement. Both models suffer from adoption of empiricism in which the principal parameters involved have little or indirect physical significance. Further development of averaged models of fingers formation is represented in [9, 26, 4]. All these models are based on the hypothesis of pressure equalization in the transverse direction to the main flow, as well as an empirical information about the displacing and displaced fluids distribution in the region of intensive viscous fingering. Let us recall that 2D numerical solution with high resolution is hard to construct, that is why 1D models play an important role in some cases. For instance, these models are very useful in calculating of fracturing when it is necessary to solve the equations of the crack opening and the fluid motion in the fracture simultaneously [1].

The aim of the present paper is to derive a hierarchy of mathematical models describing viscosity-stratified flow and spreading and growth of the mixing zone in a Hele-Shaw cell. In Section 2 we propose 2D nonlinear hyperbolic system of balance laws. In contrast with widespread model of flow displacement in a Hele-Shaw cell we apply nonlinear momentum equations and take into account compressibility of the fluid. At the same time diffusion coefficient is neglected that corresponds to the large Peclet number limit. As we show in Section 3 by numerical calculations, this model is suitable for describing of jet flow and propagation of viscous fingers in a Hele-Shaw cell. We also point out that for the process of unidirectional displacement the pressure variation in the transverse direction is small. This observation makes it possible to use long-wave approximation and construct a class of layered flows described by a system of one-dimensional evolution equations. Based on the various simplifications of the momentum equation (linearisation, lubrication theory) a hierarchy of 1D mathematical models is constructed in Section 4. The equations of a three-layer flow are studied and numerical computations of the formation of viscous fingers are performed in Section 5. We show that three-layer stationary flow is correctly described in the framework of simplified model. Nevertheless the growth rate of viscous fingers is significantly higher than it is observed experimentally. In section 6 we propose 1D kinematic-wave model of viscosity-stratified flow taking into account friction between the fluid layers. The velocity of propagation and the thickness of the viscous finger in the framework of this kinematic-wave model coincide with the corresponding calculations on the basis of the 2D equations. It gives possibility to predict the parameters of viscous fingers without time-consuming calculations. We also show that the proposed 1D model is in good agreement with the well-known Koval model.

## 2 Mathematical model

A Newtonian weakly-compressible flow displacement in a Hele-Shaw cell (the area between two parallel plates separated by a small gap of constant thickness in the direction, see Fig. 1) is described by the equations

 (ρv)t+div(ρv⊗v−P)=0,ρt+div(ρv)=0,v∣∣z=±b/2=0. (1)

Here is the fluid velocity, is the density, and is the stress tensor. This tensor can be taken in the form , where is the pressure, and is the viscosity. To describe the process of displacement involving two fluids of different constant viscosities we should take into account that viscosity depends on the concentration of solvent . This function is scaled such that it is equal to unity in the displaced fluid () and zero in the displacing one (). Following [3] we assume a monotonic relationship between the viscosity and the concentration in the form . The concentration satisfies to the transport equation (diffusion is neglected that corresponds to the large Peclet number limit)

 ct+v⋅∇c=0. (2)

We assume now that the velocity field can be represented as

 u=32(1−(2zb)2)u′(t,x,y),v=32(1−(2zb)2)v′(t,x,y),w=0. (3)

It provides the fulfilment of no-slip conditions on the cell walls . We also suppose that the functions , , and do not depend on . Let us note that in the calculation of the following terms , , and can be omitted as they are negligible compared to the derivatives with respect to . Further, averaging Eqs. (1) and (2) through the gap we obtain (primes are omitted)

 (ρu)t+(βρu2+p)x+(βρuv)y=−μu,(ρv)t+(βρuv)x+(βρv2+p)y=−μv,ρt+(uρ)x+(vρ)y=0,(cρ)t+(ucρ)x+(vcρ)y=0. (4)

Here and below denotes the viscosity of the fluid divided by the permeability (further in the text we will call it simply “viscosity”); coefficient is equal (this factor comes from integration of in the form (3) with respect to [10, 8]).

In order to close model (4) we should specify either the equation of state (barotropic fluid) or dependence (incompressible fluid). A weak compressibility of the fluid given by the equation of state provides hyperbolicity of the model. This dependence can be considered as a regularization of equations describing the flow of an incompressible fluid in a Hele-Shaw cell. On the other hand in applications the property of hyperbolicity of equations can be related to the presence of gas cavities in the porous medium as well as with the elasticity of the channel walls (for instance in PKN model [1] the pressure depends on the channel thickness ). In any case if the condition holds then results weakly depend on the choice of . Therefore for the numerical simulation of 2D flows we assume

 p(ρ)=a2ρ2/2(a2=c20/ρ0) (5)

where the constants and specify characteristic density and speed of sound in the fluid.

To find the characteristics of system (4), (5) we write it in the vector form

where is the vector of dependent variables; is the right-hand side; and are matrices. Let be the normal vector to the characteristics; is the identity matrix. Then the characteristic matrix of system (4), (5) has the form

 C=⎛⎜ ⎜ ⎜ ⎜⎝χ1+(β−1)uξ2(β−1)uξ3((β−1)(χ2−ξ1)u+p′(ρ)ξ2)ρ−10(β−1)vξ2χ1+(β−1)vξ3((β−1)(χ2−ξ1)v+p′(ρ)ξ3)ρ−10ρξ2ρξ3χ20000χ2⎞⎟ ⎟ ⎟ ⎟⎠.

Here , . A simple but cumbersome calculation yields the following expression for :

 detC(\boldmathξ)=((ξ21+2β(uξ2+vξ3)ξ1+β(uξ2+vξ3)2)−(ξ2+η2)p′(ρ))χ1χ2.

We specify the characteristic surface by the equation . Then to obtain the differential equations of the characteristics we should replace the vector in previous equation by the vector and equate to zero . As a result we obtain two families of contact characteristics

 Wt+uWx+vWy=0,Wt+βuWx+βvWy=0

 Wt+βuWx+βvWy=±√β(β−1)(uWx+vWy)2+(W2x+W2y)p′(ρ).

If the inequalities and hold this system of equations is hyperbolic. Note that in the case of , , and system (4), (5) coincides with the well-known shallow water equations.

## 3 Modelling of viscous fingering and jet flows

Below we present the results of numerical calculations of the viscous fingers and jet streams on the basis of hyperbolic model (4), (5). Originally mathematical description of the Saffman–Taylor instability was given in the framework of the Darcy’s law and the mass conservation equation [20, 22]. The inertia of the fluid may be significant for high finger velocities. In [7] simulation of viscous fingers was performed using non-linear equations of an incompressible fluid. We show that Eqs. (4), (5) taking into account the forces of inertia and compressibility of the fluid could be also used for description of this instability.

To solve differential balance laws (4), (5) numerically one can apply methods based on various modifications of Godunov’s scheme. In this work we implement the robust and stable Nessyahu–Tadmor second-order central scheme [17]. In every test we assume that on the boundaries and the impermeability condition is fulfilled. The size of the computational domain is , ; the resolution of the problem on the and axes are 300 and 150 nodes correspondingly (uniform grid). We assume that , and (for the third test ). The values of the variables are considered as dimensionless. Below we present calculations showing the possibility of modelling the evolution of perturbations caused by the Kelvin–Helmholtz and/or Saffman–Taylor instabilities on the basis of the hyperbolic model (4), (5).

### 3.1 Test 1. Jet flow

Let at the initial time the Hele-Shaw cell be occupied by a quiescent fluid having density and viscosity . Through the left central cross-section of width fluid of viscosity is injected with velocity ; through the rest part of the left boundary fluid of viscosity is entered with velocity . On the right boundary of the domain the condition of constant pressure is valid. For more intensive development of the perturbations at each time step we slightly disturb boundary conditions at . Namely, the cross section, through which fluid “2” is injected, is randomly shifted up/down from its initial position on fixed distance which is equal to the grid spacing with some positive integer factor . We call it “random shake” and take here . The function for values of the concentration is presented in Fig. 3 at . Vortices are formed at the interface of the layers due to Kelvin–Helmholtz instability. Increase of the both values and suppresses this instability.

### 3.2 Test 2. Viscosity-stratified flow

Viscosity-stratified multi-layer flow is shown in Fig. 3 at . Initially the flow region is filled by a quiescent fluid (, ). Liquid is pumped through the left boundary divided into layers of height ; velocity and viscosity for odd and even layers are , and , , correspondingly. As in the previous example “random shake” of the jets is used. The results of the calculations show that multi-layer flow without mixing is realized for a wide range of parameters (Kelvin–Helmholtz instability at the interfaces between the layers occurs if viscosity decreases more than in five times). We note that the condition on the left boundary corresponds to a class of exact solutions of equations (4) for an incompressible fluid: , , , , . The stability analysis of this class of flows was carried out in [6].

### 3.3 Test 3. Viscous fingering

The following example illustrates the formation of viscous fingers. Let the displacing phase injected at a constant velocity be referred to with index 1 and the displaced one with index 2. At fluid “1” is located in the domain ; more viscous fluid “2” — in the domain . For convenience we use the coordinate system moving with velocity and assume (with results are similar). Let us perturb the initial interface as follows: (here ). We note that the physical effect of the instability can be obtained numerically if the initial perturbation is not less than the grid resolution. At we choose piece-wise linear pressure distribution ( for and for , on the right boundary is equal to ). Initial density of the fluid is determined using formula (5). At the boundaries the impermeability condition is fulfilled.

The calculations are performed for the following parameters: , , and . In the evolution process of the flow viscous fingers are formed (Fig. 4, left). The number of fingers is determined by the initial perturbation. Displacing fluid penetrates more rapidly into displaced one (fingers are not symmetrical with respect to the initial interface). Fig. 4 (right) shows the distribution of the density at . Calculation with better resolution leads to the same result. As we can see the density changes less than in comparison with the initial one (to reduce the compressibility we should increase the speed of sound but this slows down the calculations since the time step is determined by the Courant number). Note that the density (pressure) varies slightly with respect to . This allows one to use approximate model, where the second momentum equation is replaced by .

## 4 Layered flows

We consider an incompressible () Hele-Shaw flow in a two-dimensional domain of rectangular geometry (dimensions and , respectively) governed by Eqs. (4). It is assumed that the flow is essentially parallel and the pressure gradient in the flow direction being independent of the transverse coordinate to leading-order in the variable . This regime is termed as the state of transverse flow equilibrium [26]. Such flows can be also considered in the framework of long wave approximation [14, 5]. Let us perform the following scaling in Eqs. (4)

 t→ε−1t,x→ε−1x,v→εv,μ→εμ.

Then we neglect terms of order . As a result we obtain the approximate model

 ut+βuux+βvuy+px=−μu,py=0,ux+vy=0,ct+ucx+vcy=0,v∣∣y=0=0,v∣∣y=H=0, (6)

where the pressure does not depend on the variable . We also suppose that on the boundaries and the impermeability condition is fulfilled.

Let us consider the class of viscosity-stratified flows

 u=ui(t,x),c=ci=const,y∈(yi−1,yi)

. In this case Eqs. (6) take the form

 uit+βuiuix+px=−μiui,hit+(uihi)x=0,(i=1,...,N)N∑i=1hi=H,N∑i=1uihi=Q. (7)

Here is the depth of th liquid layer of viscosity having velocity ; and is the total flow rate through the cell. Upon derivation of Eqs. (7) the kinematic condition at the layers interface is used.

Introducing new unknown variables allows to transform Eqs. (7) to the evolution system of equations

 sit+β((si/2+uN)si)x=(μN−μi)uN−μisi,hit+((si+uN)hi)x=0,(i=1,...,N−1)

where

 hN=1−N−1∑i=1hi,uN=1H(Q−N−1∑i=1sihi).

Further we assume that .

We also use the following simplified versions of governing Eqs. (7). The first one consists in the linearisation of the momentum equations:

 uit+βUuix+px=−μiui(i=1,...,N)

(here is the average velocity). The second simplification is based on the Darcy law:

 px=−μiui(i=1,...,N).

The remaining equations of system (7) do not vary.

In some cases it is convenient to use a moving coordinate system , . Then Eqs. (7) take the form (primes are omitted)

 uit+(βui+(β−1)U)uix+px=−μi(ui+U),hit+(uihi)x=0,N∑i=1hi=H,N∑i=1uihi=0. (8)

Further we show that in the framework of three-layer and two-layer regimes of flow it is possible to give an interpretation of the Saffman–Taylor instability as well as to describe the initial stage of viscous fingering.

## 5 Three-layer flow

We introduce the following notation for the layer velocities and depths

 u=u1, v=u2, w=u3; h=h1, η=h2, ζ=h3.

We also assume that , , , and .

### 5.1 Stationary solutions

Here we construct a steady-state solution of model (7) for three-layer flow. Integration of the equations of conservation of mass in system (7) allows to express the depths of the layers

 h=Q1/u,η=Q2/v,ζ=Q3/w. (9)

Here is the flow rate in the -th layer (). Due to the unit depth we obtain the velocity in the intermediate layer

 v=φ(u,w)=Q2Δ,Δ=1−Q1u−Q3w.

Eliminating pressure from the equations

 βuu′+p′=−u,βvv′+p′=−μv,βww′+p′=−w

(here the prime denotes the derivative with respect to ) reduces the problem to the solution of the autonomous system of ordinary differential equations

 dudx=(u−μφ)w−(u−w)φφw((φφu−u)w+uφφw)β,dwdx=(u−μφ)u+(u−w)(φφu−u)((φφu−u)w+uφφw)β. (10)

A fixed point of system (10) is determined from the relations :

 u∗=w∗=1+(μ−1)Q2.

Linearisation of Eqs. (10) on the solution , and computation of the eigenvalues of the corresponding matrix show that the fixed point is a stable node. The integral curves in the phase plane in the neighbourhood of the fixed point are shown in Fig. 6 for , , , and .

Fig. 6 shows a comparison of the numerical results obtained on the basis of two-dimensional hyperbolic equations (4), (5) and multilayer model (7) reduced to dynamical system (10) in the case of three-layer stationary flow. Solid white lines indicate the layers depths and obtained by solving equation (10) and using relations (9); dotted lines correspond to the fixed point (, ). Here we take the following values of layers depths , at . As before we choose , , , and . The figure shows that the solution reaches an equilibrium state for .

To carry out the calculation on the basis of 2D equations (4), (5) the following initial data are used. The flow domain in the -direction is divided into three layers of width , , and . The fluid of density at moves in these layers in the -direction with constant velocities , , and respectively. We also suppose that in the layers of width and ; in the middle layer of width we choose . The same data are taken as the boundary conditions at ; on the right boundary () the condition of constant pressure is prescribed; on the walls and the condition of impermeability is fulfilled. The resolution of the problem on the and axes are 300 and 60 nodes correspondingly (uniform grid). In order to visualize the flow of fluids with different viscosities the concentration is used. This value is presented in Fig. 6 at . More viscous fluid in the middle layer is shown in brown () and less viscous one is shown in blue ().

### 5.2 Non-stationary solutions

Let us consider a three-layer flow governing by Eqs. (7) wherein the momentum equations are replaced by linear Darcy laws . Taking into account assumptions above and notations we have

 u=w=μ/d,v=1/d,d=(1−μ)η+μ. (11)

In this case the depths of the layers and are found from the system of equations

 ht+(u(η)h)x=0,ηt+(v(η)η)x=0. (12)

It is easy to check that this system is hyperbolic and its characteristic velocities are

 λ1=ηv′(η)+v(η),λ2=u(η).

The first family of characteristics is genuinely nonlinear whereas the second one is linearly degenerate [19]. In terms of the Riemann invariants and Eqs. (12) take the form

 ηt+λ1(η)ηx=0,rt+λ2(η)rx=0.

In the case we construct a centred simple wave solution defined by the relations

 r=h0=const,λ1(η)=ξ,ξ=(x−x0)/t

(note that the ansatz leads to a constant solution). The layer depths are

 η(ξ)=11−μ(√μξ−μ),h(ξ)=(1−η)h0(μ<ξ<1μ). (13)

Formulae (13) give the solution of Eqs. (12) with discontinuous initial data

 (h,η)∣∣t=0={(0, 1),xx0. (14)

Profiles of a viscous finger and given by (13) are shown in Fig. 8 for and various values of .

Let us construct a solution of Cauchy problem (12), (14) for . At initial time the velocities and the layers depths are ,  ,  ,  ;  ,  ,  ,  . It is easy to verify that these values satisfy the Hugoniot conditions

 [(u−D)h]=0,[(v−D)η]=0

derived from Eqs. (12) as well as the stability conditions [19] if the shock front moves with average flow velocity .

It is interested to note that in the class of simple wave solutions () system (12) reduces to the naïve Koval model [4, 26]. In fact, taking into account representation (11) the second equation in (12) can be rewritten as

 ∂¯c∂t+∂∂x(M¯cM¯c+1−¯c)=0 (15)

where stands for and . It is known that the growth rate of viscous fingers in the framework of the naïve Koval model is significantly higher than it is observed experimentally.

To derive another model of a three-layer flow in a moving coordinate system we use linearisation of momentum equations in (8)

 uit+γUuix+px=−μi(ui+U),(γ=β−1).

Eliminating the pressure leads to the system of evolution equations

 s1t+γs1x=−s1,s2t+γs2x=−μs2+(1−μ)(1−w),ht+(uh)x=0,ηt+(vη)x=0 (16)

where

 u=s1+w,v=s2+w,w=−hs1−ηs2.

Let us rewrite Eqs. (16) in the form . Here is the unknown vector, is the Jacobi matrix, and is the right-hand part. The eigenvalues of matrix are

 λ1=γs1,λ2=γs2,λ3,4=2−1((1−3h)s1+(1−3η)s2±√m)

where . Conditions and provide hyperbolicity of system (16) since . Introducing the parabolic with respect to function it is easy to check validation of inequalities and (here is a minimum point of ). It means that and, consequently, characteristic velocities are real.

Further we construct the numerical solution of Eqs. (16) with initial data

 (s1,s2,h,η)∣∣t=0={(0,1−μ,0, 1),xx0.

Note that this formulation corresponds to Cauchy problem (12), (14). Calculations on the basis of model (16) are carried out using Nessyahu–Tadmor scheme [17]. The results of computations are shown in Fig. 8 at various moments of time (solid curves). As we can see with increasing of time the solution tends to the self-similar regime (dot-dashed curves obtained by using formulae (13)). Moreover tip of the viscous finger propagates with the same velocity in the framework of models (16) and (12).

## 6 Kinematic-wave model

In the previous section it is shown that three-layer flow governed by the simplified 1D model (12) (or (16)) correctly describes the well-known fact that the fluid interface is unstable if less viscous fluid displaces more viscous one and stable otherwise. However, the velocity of the viscous fingers propagation for these equations is the same as for the naïve Koval model which vastly over-predicts the rate of the fingers grow [4, 26]. A number of empirical models is proposed to reconcile this behaviour [13, 9, 24]. For example, Koval [13] postulates empirically that (15) is valid if the viscosity ratio is replaced by an effective viscosity ratio , where

 Me=(M1/4ce+1−ce)4,ce=0.22. (17)

Despite the fact that this model appears to work in practice [15] a theoretical justification is yet to be obtained.

As it can be seen from Fig. 8 and 8, in the framework of model (12) or (16), the tip of the viscous finger becomes infinitely thin with time. However, in experiments the formation of fingers of finite thickness is observed. The friction between the fluid layers is one of the possible mechanisms that prevent thinning of the viscous finger tip. Below we obtain modification of the above considered 1D models, which takes into account friction between the fluid layers. This model is new and correctly describes growth of the viscous fingers. It is proved by comparing with the Koval model (Eq. (15) where an effective viscosity ratio is used) as well as by numerical calculations on the basis of 2D Eqs. (4), (5).

### 6.1 Friction between the layers

Let us consider a two-layer flow governed by the following system of equations

 ut+βuux+px=−μ1u−κ¯μ(u−v)h−1,vt+βvvx+px=−μ2v+κ¯μ(u−v)η−1,ht+(uh)x=0,ηt+(ηv)x=0,h+η=1,uh+vη=1. (18)

Here we use notations of the previous section for the layers velocities and depths; constant is an empirical parameter, and . This system differs from Eqs. (7) (in the case ) by the presence of additional terms with factor which expresses friction on the layers interface. Note that this effect becomes significant if the thickness of one of the layers tends to zero.

Further we use simplification of the momentum equations in (18) based on the Darcy law:

 px=−μ1u−κ¯μ(u−v)h−1,px=−μ2v+κ¯μ(u−v)η−1.

Taking into account that

 η=1−h,v=(1−uh)η−1,μ1u+κ¯μ(u−v)h−1=μ2v−κ¯μ(u−v)η−1

we obtain the following kinematic-wave model

 ht+(Φ(h))x=0,Φ≡uh=(κ√M+(1−h)hM)hκ√M+(1−h)(1−(1−M)h)h, (19)

where . Typical graph of the flux is given in Fig. 10 (curve 1) for , and . Obviously is a non-convex monotonic function such that , and .

We construct self-similar solution of Eq. (19) with initial data for and for . On the interval function has two points of inflection. Thus it is necessary to construct convex hull of the flux . For this we draw tangents to from the origin and from the point (lines 2 ans 3 in Fig. 10). Let and be the tangency points (values and are obtained by solving the equations and ). According to [19] solution of the Cauchy problem takes the form

 ξ=⎧⎪⎨⎪⎩Φ′(h1),h∈[0,h1)Φ′(h),h∈[h1,h2]Φ′(h2),h∈(h2,1]

where is the self-similar variable. This means that the solution is presented by centred rarefaction wave, which is bounded by two “sonic” shocks. As we show below the numerical results of the growth of viscous fingers in the framework of Eqs. (4), (5) are in good agreement with the self-similar solution of Eq. (19) for a suitable choice of the parameter .

### 6.2 Comparing with the Koval model and numerical results

Let us verify kinematic-wave model (19) by comparing with the well-approved Koval model (15), where instead of effective viscosity ratio given by formula (17) is used. We choose the following parameters: , . It means that the viscosity ratio and according to (17) . Self-similar solution of the Koval model for these parameters is shown in Fig. 10 (solid curve). Corresponding solution of the scalar conservation law (19) (with ) is presented in Fig. 10 by dashed curve. As we can see in the framework of these models the growth rates of viscous fingers are similar. Although the parameter is a function of , the values of are weakly vary for . Therefore, for the moderate viscosity ratios (less than 10) we can assume that . Moreover, the proposed model (19) describes better the structure of the finger (its thickness near the tip) than the Koval model (15).

Further we compare the results of numerical simulations of the formation of viscous fingers on the basis of 2D hyperbolic system of equations (4), (5) with the results obtained by using kinematic-wave model (19). The calculations are performed in the Hele-Shaw cell with sizes , in the coordinate system moving with the average flow velocity with respect to axis. The viscosities of the fluids are equal to and . At the initial time the interface is perturbed as follows , where . For discretization with respect to and we use 400 and 50 nodes respectively (calculation with finer resolution leads to the same result). On the left and right boundaries of the computational domain the reflection conditions are imposed. In 2D model (4), (5) the following dependence for viscosity is used: . The calculations are carried out for , and the sound velocity . In this case the change in density is not more than 0.15%. At the same time the condition is fulfilled with high accuracy. It means that the above proposed 1D model is suitable for describing of such flow. Given above perturbation of the interface leads to the formation of a single finger which is symmetric with respect to the line . The results of the concentration calculations using model (4), (5) at time and are shown in Fig. 11. The distribution of is presented in two colours (blue for and brown for ).

Self-similar solution of kinematic-wave model (19) with is shown in Fig. 11 at and by white solid lines respectively. The curves and for the model are given for a correct comparison with 2D calculations. Here self-similar variable is replaced by . It corresponds to a transition in a moving coordinate system. The figure shows the propagation velocity and thickness of the viscous finger obtained by the kinematic-wave model agree well with the two-dimensional calculations.

## 7 Conclusion

We derive nonlinear hyperbolic system of equations (4), (5) describing the flow of slightly compressible multicomponent fluid of different viscosity in a Hele-Shaw cell. On the basis of these equations simulation of jet flow and development of viscous fingers during the displacement process are performed. Calculations show that the proposed model reproduces the characteristic features of the flow associated with the development of Kelvin–Helmholtz and Saffman–Taylor instabilities (see Fig. 3, 3 and 4). In the case of preferential flow in the -direction the pressure varies slightly in the -direction that allows to apply model of long-wave approximation (6) and to consider the class of layered flows described by Eqs. (7). Using various simplifications of the system (linearisation of the momentum equations, application of the Darcy’s law) we construct a hierarchy of mathematical models of viscosity-stratified flow in a Hele-Shaw cell. These 1D models are suitable for description of the main features of the two-dimensional flow. Stationary solutions of Eqs. (7) obtained for the three-layer flow are in good agreement with the calculations of the flow on the basis of 2D model (4), (5) (see Fig. 6). In the framework of the three-layer flow (systems (12) and (16)) the interpretation of the Saffman–Taylor instability is given. Solutions of Eqs. (12) and (16) illustrating the formation of viscous fingers are constructed (Fig. 8 and 8). However, these models fail to predict the correct growth rate of the fingers.

We propose modification of the layered flow model in order to agree with this behaviour. The model is obtained by including friction between the fluid layers (18). This provides a non-zero thickness of the fingertip and under some assumptions allows one to describe the evolution of viscous fingers on the basis of the scalar equation with non-convex flux (19). Although the proposed equation (19) involves empirical parameter, this model reveals the physical mechanism to ensure correct propagation velocity and structure of viscous fingers. As it can be seen from Fig. 11 the velocity of propagation and the thickness of the fingers in proposed model (19) (if the empirical parameter is appropriately specified) are in fairly good agreement with calculations based on two-dimensional equations (4), (5). Comparison with the Koval model (see Fig. 10) confirms the correctness of the results.

## Acknowledgements

This work was supported by the Russian Science Foundation (grant No. 15-11-20013).

### References

1. Adachi, J., Siebrits, E., Peirce, A. & Desroches, J. (2007) Computer simulation of hydraulic fractures. Int. J. Rock Mech. Min. Sci. 44, 739–757.
2. Al-Housseiny, T. T., Tsai, P. A. & Stone, H. A. (2012) Control of interfacial instabilities using flow geometry. Nat. Phys. 8, 747–750.
3. Azaiez, J. & Singh, B. (2002) Stability of miscible displacements of shear thinning fluids in a Hele-Shaw cell. Phys. Fluids 14, 1559–1571.
4. Booth, R. J. S. (2010) On the growth of the mixing zone in miscible viscous fingering. J. Fluid Mech. 655, 527-539.
5. Chesnokov, A. A. & Liapidevskii, V. Yu. (2011) Shallow water equations for shear flows. Notes Numer. Fluid Mech. Multidisciplinary Design. 115, 165–179.
6. Chesnokov, A. A. & Stepanova, I. V. (2015) Stability analysis of shear flows in a Hele-Shaw cell. Appl. Math. Comput. 265, 320–328.
7. Chevalier, C., Amar, M., Bonn, D., & Linder, A. (2006) Inertial effects on Saffman — Taylor viscous fingering. J. Fluid Mech. 552, 83–97.
8. Dias, E. O. & Miranda, J. A. (2011) Influence of inertia on viscous fingering patterns: rectangular and radial flows. Phys. Rev. E. 83, 066312.
9. Fayers, F. J. (1988) An approximate model with physically interpretable parameters for representing viscous fingering. SPE Reservoir Engng 5, 551–558.
10. Gondret, P. & Rabaud, M. (1997) Shear instability of two-fluid parallel flow in a Hele-Shaw cell Phys. Fluids 9, 3267–3274.
11. Govindarajan, R. & Sahu K. C. (2014) Instabilities in viscosity-stratified flow. Annu. Rev. Fluid Mech. 46, 331–353.
12. Homsy, G. M. (1987) Viscous fingering in porous media. Annu. Rev. Fluid Mech. 19, 271–311.
13. Koval, E. J. (1963) A method for predicting the performance of unstable miscible displacement in heterogenous media. Soc. Petrol. Eng. J. 450, 145–154.
14. Liapidevskii, V. Yu. & Teshukov, V. M. (2000) Mathematical Models of Propagation of Long Waves in a Non-Homogeneous Fluid. Novosibirsk. Siberian Division of the Russian Academy of Sciences. (in Russian).
15. Malhotra, S., Sharma, M. M. & Lehman, E. R. (2015) Experimental study of the growth of mixing zone in miscible viscous fingering. Phys. Fluids 27, 014105.
16. Manickam, O. & Homsy, G. M. (1991) Fingering instabilities in a vertical displacement flows in porous media. J. Fluid Mech. 288, 75–102.
17. Nessyahu, H. & Tadmor, E. (1990) Non-oscillatory central differencing schemes for hyperbolic conservation laws. J. Comp. Phys. 87, 408–463.
18. Pihler-Puzovic, D., Perillat, R., Russell, M., Juel, A. & Heil, M. (2013) Modelling the suppression of viscous fingering in elastic-walled Hele-Shaw cells. J. Fluid Mech. 731, 162–183.
19. Rozhdestvenskij, B. L.& Yanenko, N. N. (1983) Systems of quasilinear equations and their applications to gas dynamics. Transl. Math. Monogr., 55, Amer. Math. Soc., Providence.
20. Saffman, P. G. & Taylor, G. (1958) The penetration of a fluid into a porous medium or a Hele-Shaw cell containing a more viscous liquid. Proc. Roy. Soc. A 245, 312–329.
21. Smirnov, N. N., Nikitin, V. F., Maximenko, A., Thiercelin, M. & Legros, J. C. (2005) Instability and mixing flux in frontal displacement of viscous fluids from porous media. Phys. Fluids 17, 084102.
22. Tan, C. T. & Homsy, G. M. (1986) Stability of miscible displacements in porous media: rectilinear flow. Phys. Fluids 29, 3549–3556.
23. Tanveer, S. (2000) Surprises in viscous fingering. J. Fluid Mech. 409, 273–308.
24. Todd, M. R. & Longstaff, W. J. (1972) The development, testing, and application of a numerical simulator for predicting miscible flood performance. J. Petrol. Tech. 24, 874–882.
25. Yang, Z. M., Yortsos, Y. C. & Salin, D. (2002) Asymptotic regimes in unstable miscible displacements in random porous media. Adv. Water Res. 25, 885–898.
26. Yortsos, Y.C. & Salin, D. (2006) On the selection principle for viscous fingering in porous media. J. Fluid Mech. 557, 225–236.
27. Yuan, Q. & Azaiez, J. (2015) Inertial effects of miscible viscous fingering in a Hele-Shaw cell. Fluid Dyn. Res. 47, 015506.
28. Zimmerman, W. B. & Homsy, G. M. (1992) Viscous fingering in miscible displacements: unification of effects of viscosity contrast, anisotropic dispersion, and velocity dependence of dispersion on nonlinear finger propagation. Phys. Fluids A 4, 2348–2359.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters