Three-dimensional aspects of fluid flows in channels. I. Meniscus and Thin Film regimes

We study the forced displacement of a fluid-fluid interface in a three-dimensional channel formed by two parallel solid plates. Using a Lattice-Boltzmann method, we study situations in which a slip velocity arises from diffusion effects near the contact line. The difference between the slip and channel velocities determines whether the interface advances as a meniscus or a thin film of fluid is left adhered to the plates. We find that this effect is controlled by the capillary and Peclet numbers. We estimate the crossover from a meniscus to a thin film and find good agreement with numerical results. The penetration regime is examined in the steady state. We find that the occupation fraction of the advancing finger relative to the channel thickness is controlled by the capillary number and the viscosity contrast between the fluids. For high viscosity contrast, Lattice-Boltzmann results agree with previous results. For zero viscosity contrast, we observe remarkably narrow fingers. The shape of the finger is found to be universal.


I. INTRODUCTION
Advancing fronts in fluid systems involve the motion of a fluid-fluid interface, a surface that lives in a threedimensional world, and which is often constrained by a solid boundary. A typical example is that of an interface moving in a channel [1,2,3].
Examples of advancing fronts in channels are imbibition, a process in which a wetting fluid invades the channel due to an uncompensated capillary pressure, and the viscous fingering process [2,3,4], where a low-viscosity(or high-density) fluid penetrates a high-viscosity(or lowdensity) one.
The problem of viscous fingering in a channel has been widely studied in the framework of Hele-Shaw theory. A Hele-Shaw cell is the two-dimensional limiting case of a very thin channel, where the equations of motion are averaged over the channel thickness. This reduces the interface to a line, the leading interface, that lives in the plane of the cell. As a consequence, the approximation discards any effects arising from the full three-dimensional structure of the interface.
Nonetheless, penetration in the gap of a Hele-Shaw cell is a fundamental three-dimensional effect that has important repercussions in the viscous fingering problem. As theoretical studies have pointed out [5], a thin film of viscous fluid left adhered to the cell plates as the front advances modifies the capillary pressure at the leading interface, thus altering the front morphology. This has been confirmed in experiments of steady viscous fingers [6], where the presence of a thin film led to fingers not predicted by the two-dimensional theory. In * Electronic address: rodrigo@ecm.ub.es the following paper, we will address the role of the thin film in viscous fingers.
A thin wetting film is not the only consequence of a three-dimensional interfacial structure. In the context of liquid films spreading over dry substrates [7,8], where a two-dimensional approximation is typically applied, three-dimensional effects are also important. For instance, the stability of a spreading front depends on the wetting properties of the fluid. Experimentally, it has been observed [9] that a crossover from unstable to stable fronts occurs when the dynamic contact angle exceeds π/2, a situation which renders the velocity field within the film three-dimensional.
The problem of a moving interface in a threedimensional channel must take into account a dynamic contact line, the intersection point between the fluidfluid interface and the channel walls. In classic fluid mechanics, a moving contact line violates the usual nostick boundary condition, leading to a divergent viscous dissipation [10]. Hence, contact line dynamics must consider a regularizing mechanism of the viscous dissipation singularity. A slip velocity in the vicinity of a driven contact line arises naturally in diffuse interface models of binary fluids [11], regularizing the singularity. These models consist of the usual Navier-Stokes equations coupled to a convection-diffusion equation of an order parameter [12]. Diffuse interface effects enter the force balance equations in the shape of order parameter gradients that play the role of a Young force. As a result, the contact line slips over the solid boundary [11,13]. Away from the contact line, order parameter gradients vanish and the stick boundary condition is recovered. The size of the diffusion region, l D , is then a measure of how strong slip is for a given system and is clearly an important parameter. This size was estimated by Briant and Yeomans [14], who characterize l D for the case of an interface subjected to shearing walls. They focused on the dependence of l D (L in their notation) on the model parameters, finding a scaling relation that was verified numerically.
Important implications arising from a relatively large or small slip velocity compared to the leading interface velocity in forced fronts can be foreseen. Whenever both velocities are comparable, the interface should maintain a meniscus shape. Conversely, as the slip velocity becomes small compared to the channel velocity the interface shape should develop as a finger, leaving a thin film of fluid adhered to the walls of the channel.
In this paper we study the penetration process across the channel thickness. We study the motion of the full three-dimensional interface between two viscous fluids when it is subjected to a gravitational body force. We treat the case of a strictly flat leading interface, focusing only in the details that pertain to the channel thickness. We work with symmetric fluids as well as with fluids of different densities or viscosities.
We focus on two principal matters. We first describe how the contact line and leading interface velocities are related, and propose the mechanisms that determine the velocity ratio. We find that the velocity ratio is controlled by the force balance at the interface and by diffusion effects localized at the contact line.
Secondly, we study the thin film that forms inevitably in the case of small slip. In that case the front decouples from the contact line leading to the growth of a finger, even when the interface is linearly stable to the Rayleigh-Taylor instability. We find that the fraction of occupation of the thin film relative to the channel thickness is a function of the capillary number and the viscosity contrast between the fluids. The high viscosity contrast case is validated by comparing our results to the numerical work of Halpern and Gaver [15] which is consistent with previous results of Taylor [16]. For fluids with zero viscosity contrast, it turns out that the finger width has much lower values than for the high viscosity contrast case at fixed capillary number.
The morphology of our fingers is very much alike to the Saffman-Taylor finger shape, a prediction of the Hele-Shaw theory. Nevertheless, it is important to stress that although the case of a flat leading interface is twodimensional, the equations of motion are not equivalent to those of the Saffman-Taylor problem. Therefore, penetration in the channel thickness cannot be attributed to the Saffman-Taylor instability. Likewise, the selection rule of the steady state, i.e., the actual dependence of the thin film thickness with the front velocity, cannot be mapped to the theoretical predictions of the viscous fingering theory.
We will address these matters by means of numerical simulations of the mesoscopic equations of the system. To do so, we take advantage of a powerful integration algorithm in fluid dynamics, the Lattice-Boltzmann scheme for binary fluids.
The paper is organized as follows. In Sec. II we present the equations that govern the system in the meso-scopic regime. In Sec. III we briefly present the Lattice-Boltzmann algorithm for binary fluid flows. Sec. IV A is dedicated to simulation results of the forced interface, from which two steady state regimes are found; a nonpenetrating regime, in which the interface advances as a meniscus, and a penetrating one, in which a single finger emerges and achieves steady state. In Sec. IV B we present a scaling argument of the equations of motion that leads to an estimate of the ratio between the slip and front velocities. Such argument explains the crossover from one regime to the other. In Sec. IV C we extend our results to fluids of different densities or viscosities. Sec. IV D is devoted to the steady state finger. Finally, in Sec. V we present the conclusions of this work.

II. GOVERNING EQUATIONS
We consider a channel formed by two solid plates parallel to the xy plane, each of length L and infinite width, located at positions z = 0 and z = b. Initially, two fluids fill the channel and are separated by a flat interface perpendicular to the solid walls, as shown in Fig. 1. The equilibrium contact angle is hence θ E = π/2. Contact lines are located at z = 0 and z = b, while the leading interface is located at z = b/2.
To circumvent the complications of the sharp interface formulation, we introduce a mesoscopic variable, φ( r); and order parameter which is constant in the bulk of each fluid and varies smoothly across a diffuse interfacial region. Within this approximation, the equilibrium state of the system is described by a Helmholtz free energy [12] The first term in the integrand is a volume contribution, given by V (φ, ρ) = Aφ 2 /2 + Bφ 4 /4 + ρ/3 ln ρ. The ρ dependent term corresponds to an ideal gas contribution, while the φ dependent terms allow for the coexistence of two phases. The presence of an interface is accounted for by the last term in the integrand, which penalizes spatial variations of the order parameter by a factor κ. Minimization of F leads to the chemical potential, and total pressure tensor [17], whereδ is the diagonal matrix. The pressure tensor has an ideal contribution given byP = ρ 3δ , and an order parameter contribution. In equilibrium, the order parameter profile for the flat interface (sketched in Fig. 1) is φ * (x, z) = −φ eq tanh(x/ξ), where φ eq = (−A/B) 1/2 is the bulk equilibrium value of the order parameter and ξ = (−κ/2A) 1/2 is the length scale of the interfacial region; this profile leads to the difference between equilibrium values ∆φ = 2φ eq and the energy per unit area of the interface, σ = (−8κA 3 /9B 2 ) 1/2 . Since the interface is diffuse, a choice for the nominal interface position has to be made. We choose the level surface φ = 0.
The divergence of the pressure tensor yields the force per unit volume that acts on the fluid: − ∇P − φ ∇µ. The first term is the pressure gradient, while the second arises from order parameter inhomogeneities. Consequently, the Navier-Stokes equations are [12], where v is the fluid velocity, η is the fluid viscosity and g is the acceleration of gravity. The dynamics of the order parameter are described by a convection-diffusion equation, where M is a mobility. For small deviations from the equilibrium configuration, an expansion of the chemical potential in powers of φ − φ * yields a first order diffusion coefficient D = M (A + Bφ 2 eq ), so the relative importance of the advective and diffusive terms can be estimated through a Péclet number, P e = | v · ∇φ|/|D∇ 2 φ|.
The system can be represented as a sheet of fluid in the xz plane with periodic boundary conditions applied in the y direction. This is equivalent to a channel of infinite width in the y direction with a flat leading interface. Stick boundary conditions are imposed at the walls, v(x, z = 0) = v(x, z = b) = 0, while no flow boundary conditions are imposed for the order parameter, φ v(x, z = 0) = φ v(x, z = b) = 0. At both ends of the channel the flow is homogeneous. This is ensured by setting ∂ Contact line dynamics arise from the diffuse nature of the interface, which allows for slip in the interfacial region by a diffusive mechanism. The size over which slip takes place, l D , is a function of the fluid properties and has been estimated by Briant and Yeomans [14], who have given a scaling relation, l D ∼ (ηξ 2 M/∆φ 2 ) 1/4 .

III. LATTICE BOLTZMANN METHOD
We solve numerically Eqs. (1) and (2) by means of the Lattice-Boltzmann algorithm presented in Ref. [17]. The dynamics are introduced by discretized Boltzmann equations of two distribution functions, and In these equations, f i and g i are distribution functions, where the index i counts over the model velocity set. Space is discretized as a cubic lattice where nodes are joined by velocity vectors, c i . Space and time units in Eqs. (3) and (4) are set to unity. Likewise, the density of the fluids is set to one. We use the D3Q15 velocity set, which consists of fifteen velocity vectors: six of magnitude 1 that correspond to nearest neighbors, eight of magnitude √ 3 that correspond to third-nearest neighbors and one of zero magnitude that accounts for rest particles. In the D3Q15 model the speed of sound is c s = 1/ √ 3. In Eqs. (3) and (4), distribution functions are first relaxed to equilibrium values, represented by f eq i and g eq i , with relaxation timescales τ f and τ g . The term F f i is related to the external forcing. Following the collision stage, distribution functions are propagated to neighboring sites. Hydrodynamic variables are defined through moments of the f i and g i . The local density and order parameter are given by i f i = ρ and i g i = φ. The fluid momentum and order parameter current, are defined as Local conservation of mass and momentum is enforced through the con- In equilibrium, the pressure tensor and chemical potential are defined as The equilibrium distribution functions and the forcing term are written as expansions in powers of v [18], i.e., Here, ν stands for the three possible magnitudes of the c i set. Coefficient values are ω 0 = 2/9, ω 1 = 1/9 and Eqs. (1) and (2) can be recovered as a Chapman-Enskog expansion of Eqs. (3) and (4) [18]. The Lattice-Boltzmann scheme maps to the hydrodynamic model through the relaxation timescales, i.e., η = (2τ f − 1)/6 and M = (τ g − 1/2)M, and through the body force f = ρ g.
Solid boundaries in the Lattice-Boltzmann method are implemented by means of the well known bounce-back rules [18,19]. In the lattice nodes that touch the solid, the propagation scheme is modified so the distribution functions are reflected to the fluid rather than absorbed by the solid. As a consequence, a stick condition for the velocity is recovered approximately halfway from the fluid node to the solid node.

IV. RESULTS
We study the process of penetration across the channel thickness in the presence of a dynamic contact line. As we have explained above, fingering is expected whenever the slip velocity is small compared to the leading interface velocity. In our model, slip is controlled by diffusion in the vicinity of the contact line. To measure the importance of diffusivity we use a typical definition of the Péclet number, P e = U b/D, where U is the velocity of the leading interface. The other relevant control parameter is the capillary number, which follows from the ratio between viscous and capillary forces, Ca = ηU/σ. We focus on flows governed by viscous and capillary forces. To enforce this situation we neglect the convective term in Eq. (1). To assure that we work on the low Mach number regime, the fluid velocity is restricted to U ≤ 0.01. For the case of small slip, we expect a thin film regime typical of experiments. We characterize this regime in terms of the finger width, viscosity contrast and capillary number. We compare our results with other studies from the literature.

A. Effect of diffusivity, surface tension and viscosity
We first consider two fluids with equal viscosities and densities. The size of the interface is set to ξ = 0.57, which has been previously verified to give sufficiently accurate results for the variation of φ and its spatial derivatives across the interface [17]. Starting from a flat interface configuration, we perform a set of five runs at fixed forcing, viscosity and surface tension. For each run we choose a different diffusion coefficient, which we fix through the mobility. In terms of dimensionless numbers this corresponds to fix Ca and vary P e. Parameter val- In Fig. 2 we show a time sequence of the interface position for each run. In our simulations, v y = 0 so a flat leading interface is located at z = b/2. Sequences (a) and (b) correspond to runs with the highest diffusion coefficients. In both cases a steady meniscus is clearly observed. It is also appreciable that the meniscus in sequence (a), corresponding to the highest diffusivity, is less curved than the meniscus in sequence (b). The next three sequences, (c), (d) and (e), show an abrupt change in the interface configuration. Instead of a meniscus, we observe a penetrating structure that emerges from the center of the channel leaving a thin film of fluid adhered to the solid plates. The finger width in runs (c)-(e) is approximately 17 lattice spacings. For the size of the interface used, the order parameter saturates to its equilibrium value at the solid surface. Nonetheless, to rule out any effects associated to the size of the interface, we have verified that the finger width (relative to the channel thickness) does not depend on b, as we will see bellow.
All runs achieve a steady state in which the velocity of the leading interface is constant. This velocity is the same for runs (a) and (b) and due to mass conservation is slightly larger (a few percent) for runs (c), (d) and (e). The capillary number is not affected much by this effect, and we will take it as constant. The relevant effect is associated to the variation of diffusivity. The velocity of the contact line increases with increasing diffusivity, as can be deduced from the contact line position in sequences (c), (d) and (e). Nevertheless, the velocity of the leading interface and the width of the penetrating finger are the same for all three runs. This is a direct confirmation of the fact that contact line dynamics are decoupled from leading interface dynamics in the presence of a thin film, as proposed by Park and Homsy in Ref. [5].
It is clear from these runs that the crossover for penetration is set by the difference between the leading interface velocity, U , and the slip velocity at the contact line, v s . For a meniscus, v s = U , while penetration occurs whenever v s < U . As v s depends on the strength of diffusivity, we can draw as a conclusion that penetration can be achieved by increasing P e.
We now explore the effect of capillarity on the dynamics of the interface. To do so we force the interface at fixed velocity, diffusivity and surface tension(resp. viscosity) while we vary the viscosity (resp. surface tension). As a consequence, P e is fixed while Ca is varied.
Results are summarized in Table I. The first column shows parameters for runs in which the viscosity is varied. We observe that penetration occurs as η increases. The second column in Table I shows results for varying surface tension. We observe that penetration occurs as σ is decreased. We can conclude that capillarity plays a similar role as diffusivity, as penetration occurs for low values of Ca.

B. The Onset of Penetration
Our results suggest that the crossover from the meniscus regime to the thin film regime is controlled at least by two mechanisms. On the one hand, viscous stresses deform the interface. As a result surface tension tends to restore the interface shape to its equilibrium value. On the other hand, advection causes order parameter gradients. As a consequence, diffusivity generates a slip velocity at the contact line. In this section we will see that the balance between these mechanisms is controlled by P e and Ca.
Let us write the force balance per unit volume of fluid in the frame of reference of the interface. We introduce orthogonal curvilinear coordinates, s, the arclength along the curve φ = 0, and u, the normal distance to a point on this curve. In terms of these coordinates the normal component of Eq. (1) (in absence of inertial terms) is where the subscript n stands for the normal component and the subscript u denotes differentiation with respect to u. The force per unit area acting on the interfacial region is obtained by integrating (5) across the interface: where the term σ(κ D σ − κ E σ ) arises from the integration of the chemical potential term [12], with κ D σ and κ E σ being the dynamic and equilibrium curvatures, which are positive for a bump protruding in the x direction. We have assumed that neither of the last two terms in the right hand side vary appreciably across the interface. Eq. (6) should be interpreted as the usual Gibbs-Thomson condition plus a dynamic term proportional to ξ, which vanishes either in equilibrium or in the sharp interface limit.
We will now examine Eq. (6) in the vicinity of the contact line. The mesoscopic nature of the interface gives rise to a finite size region where diffusion is important. This results in a slip velocity, v s , for the contact line. We now reproduce the scaling argument presented in Ref. [14] to obtain the diffusion size, l D , and consequently v s . We will subsequently compare the slip velocity to the leading interface velocity in terms of Ca and P e, which are parameters that can be linked to experiments.
The slip velocity and the size of the diffusion region fix the magnitude of viscous dissipation in Eq. (6), Since in the contact line region the time variation of the order parameter is ∂φ/∂t ≃ v s ∆φ/ξ, the order parameter variation obeys Using Eq. (8) to eliminate l D from Eq. (7) we get The last term in this expression is order ξ, while the rest of terms are of order ξ 0 . The term in the left hand side is the excess pressure drop caused by the curvature difference, which is small in our simulations. Neglecting both the pressure gradient and the body force we extract the following scaling law for the slip velocity: The interface curvature is a consequence of the underlying velocity profile, which is set by the thickness of the channel. Therefore, the curvature difference scales as (κ D σ −κ E σ ) ∼ a/b, with a being a typical amplitude. Using this expression and measuring v s in units of the leading interface velocity we arrive at the following expression: This indicates that both P e and Ca control how the slip velocity compares to the leading interface velocity. For a meniscus v s = U , so we arrive at the following condition: To check the validity of the prediction we perform several runs of forced interfaces varying simulation parameters in a wide range(see Tables III and IV). We cover up to four decades in the P e and Ca until numerical stability issues of the code begin to show up. In Fig. 3 we show a plot of our results in the Ca −1 P e plane. Data shown in the figure sketches two regions; at high P e and Ca values the thin film regime is observed, whereas the meniscus corresponds to low values of both parameters. We also show our prediction, for which we find a fitting value for a, namely a ≃ 0.3.
Typical experiments with molecular liquids correspond to high, O(10 2 ) − O(10 3 ), values of the Péclet number [6,20], thus falling in the thin film region sketched in the diagram. However, menisci are expected in systems with a diffuse interface, such as colloid-polymer mixtures [21]. In terms of experimental parameters, the diffusion length scales as l D ∼ (ηξ 2 D/σ(κ D σ −κ E σ )) 1/4 . In colloid-polymer mixtures the ratio (ξ 2 /σ) 1/4 is about 10 2 times larger than for molecular liquids. As a consequence, in such systems menisci should be observable at relatively high capillary numbers. For molecular liquids, this effect can be achieved by decreasing the system size, for instance, in microfluidic devices, where the typical size of the channel is a few microns [22]. For such small sizes, the Péclet number is O(1), about 10 3 times smaller than for traditional channels. Hence, the transition from menisci to thin films would be observable in the regime of relatively high capillary number, say Ca O(10 −1 ) − O(10 0 ).

C. Asymmetric Fluids
We have shown that the ratio between the leading interface and contact line velocities is controlled by the interplay between the Péclet and capillary numbers. We expect that this fact holds for fluids of either different densities or viscosities. A forced interface between asymmetric fluids can be destabilized by virtue of the Rayleigh-Taylor instability, when the more dense fluid displaces the less dense one. We explore situations for which the instability is absent. To ensure this, we keep the channel thickness below the first unstable wavelength, given by l c = 2π(σ/∆ρg) 1/2 [23].
The forcing is set according to where the local viscosity is set according to the mixing rule η(φ) = (η 2 + η 1 )(1 − cφ/φ eq )/2, characterized by the viscosity contrast c = (η 2 − η 1 )/(η 2 + η 1 ), and U exp is the maximum expected velocity for a Poiseuille profile. The φ dependent part is set as A(φ) = 1 if c = 0 and A = (φ + φ eq )/∆φ otherwise. In experiments the typical situation is that an effectively inviscid fluid displaces a viscous one, which corresponds to c → 1. We approach this situation by setting c = 0.9. Following the general convention in the literature, here we define the capillary number as Ca = η 2 U/σ. For all cases, we consider that both fluids have the same diffusion coefficient.
We have performed a set of runs in which we vary P e at fixed Ca for fluids with finite density or viscosity contrast. The details of the runs are summarized in Table II. For each case, both menisci and fingers can be obtained depending on the value of the Péclet number. As expected, penetration occurs for sufficiently high P e. We can conclude that the appearance of a thin film is independent of the Rayleigh-Taylor instability. In Fig. 3 we plot results of this section in the P eCa −1 plane. For fluids of different densities this value is consistent with the symmetric estimate of a ≃ 0.3. For fluids of different viscosities the crossover occurs at a ≃ 0.5. Anyhow, the qualitative behavior remains independent of the degree of asymmetry between the fluids.

D. Steady state finger in the channel thickness direction
We now turn our attention to the steady state finger that appears for high values of the product CaP e, which is the usual situation in most experiments. As we have shown in Sec. IV A, diffusion only affects the motion of  In Fig. 4 we plot λ b as a function of Ca. We find that λ b depends on the viscosity contrast, as the c = 0 points fall in a clearly different curve than the c = 0.9 and c = 0.95 points.
On the contrary, the density contrast does not play a relevant role. Points belonging to the c = 0 curve were obtained using five different gap sizes; b = 147, b = 23, b = 35, b = 51 and b = 95, of which the b = 35 set was done with fluids of different densities. Results show no difference between zero or finite density difference. In fact, the gap size for the latter was large enough for the Rayleigh-Taylor instability to be present. This means that the finger can develop as a consequence of low diffusion or by virtue of the Rayleigh-Taylor instability. Still, the steady state remains insensitive to the mechanism that leads to penetration and is selected by Ca and c.
Previous analytic predictions correspond to the low Ca regime at c = 1 and were carried out first by Bretherton [24], who found that the finger width decays as λ b → 1 − 1.337Ca 2/3 , as Ca → 0. An extension was done by Taylor [16], up to Ca < 0.09, for which he reported a decaying exponent of one-half. Numerically, Reinelt and Saffman [25] solved the Stokes equations using a finite difference algorithm, and considered values up to Ca < 2, which match the one-half exponent of Taylor at small Ca. Halpern and Gaver [15] extended the prediction beyond Ca = 2 by means of a boundary element analysis of the Stokes equations. Their results can be fitted to an exponential law λ b = 1 − 0.417 1 − exp(−1.69Ca 0.5025 ) (shown in Fig. 4), which reproduces Reinelt and Saffman results and matches the power law prediction of Taylor for low Ca. For large Ca, this law saturates to a limiting value of λ b = 0.583. Previous Lattice-Boltzmann studies have also addressed this problem. Kang et al [26] studied the range 0.2 ≤ Ca < 2, obtaining good agreement with Halpern and Gaver results. Langaas and Yeomans [27] considered the range 0.079 ≤ Ca ≤ 4.6 and were able to reproduce the results of Halpern and Gaver for Ca < 2. For Ca > 2 they obtained smaller finger widths than those of Halpern and Gaver.
Our results cover up to five decades in the capillary number, 10 −2 ≤ Ca ≤ 10 2 [28]. They match Halpern and Gaver prediction as c → 1. Already at c = 0.90, we reproduce accurately the finger width saturation value, for which we find λ b = 0.573 ± 0.022. For small Ca, the error increases for the c = 0.9 runs. We improve this situation by increasing the viscosity contrast, as can be appreciated in Fig. 4. At Ca ≃ 0.09 the error for c = 0.9 is 4%,while for c = 0.95 it reduces to 2%. For Ca = 0.008 the error is 2% at c = 0.95. This agreement shows that the Lattice-Boltzmann approach gives accurate results for a wide range of Ca, improving previous results [27].
As can be seen from Fig. 4, fingers with zero viscosity contrast, a case that has not been studied previously, are much narrow than fingers with c = 0.9 or c = 0.95. The dependence of the finger width with Ca has a power law behavior for 0.1 ≤ Ca ≤ 1, with an exponent m = 0.29 ± 0.02. For Ca O(10), the finger width saturates to a notably small value, λ b ≃ 0.386 ± 0.014, which remains an open question.
We now focus on the shape of the steady finger. Finger profiles shown in Fig. 2 are strongly reminiscent of the single finger solution of the Saffman-Taylor problem [4], in which fingering occurs in the xy plane of a Hele-Shaw cell as a consequence of viscosity or density asymmetries between the fluids. Our problem is fundamentally different because penetration in the channel thickness can occur for linearly stable interfaces in the context of hy-drodynamic stability, e.g., by virtue of low diffusivity at the contact line. Moreover, even in the case where the interface is linearly unstable, it is due to a Rayleigh-Taylor instability and not through the Saffman-Taylor one.
Still, we compare our finger profiles with the Saffman-Taylor ones. To do so, we recall the semi-empiric shape found by Pitts [29], which for our geometry reads where x ′ and z ′ measure the distance from the finger tip. For both c = 0 and c = 0.9, Eq. (10) is a good approximation to our profiles at low Ca. This agreement is lost gradually as Ca increases. In Fig. 5

V. CONCLUSIONS
We have studied the forced motion of a fluid-fluid interface in a three-dimensional channel by means of a mesoscopic model that takes into account contact line dynamics.
Our results describe two possible scenarios regarding interface dynamics. A meniscus regime is found whenever the contact line velocity is comparable to the leading interface velocity. Conversely, when the contact line velocity is smaller than the leading interface velocity the meniscus configuration is lost, leading to penetration of one fluid into the other in a fingering fashion. A thin film of displaced fluid is hence left adhered to the plates of the channel.
The crossover from the meniscus to the thin film regime is controlled by the competition between surface and viscous stresses, as well as by the competition of diffusive and advective timescales, on top of the usual hydrodynamic instabilities. These mechanisms can be accounted for through simple scaling arguments. We find a prediction for the crossover in terms of the capillary and Péclet numbers which describes accurately our numerical results. Menisci are found for low capillary and Péclet numbers, when surface tension and diffusion dominate over viscous stresses and advection respectively. An example of a system where diffusion is important is that of colloid-polymer mixtures. For such systems, the relatively large size of the interface together with low surface tensions leads to large diffusion regions near de contact line. For instance, in Ref. [21] the size of the interface is typically ξ ≃ 10 µm while the surface tension is σ ≃ 1µN. For molecular liquids the size of the interface is of the order of nanometers, while σ ≃ 10 mN/m. As explained in Sec.IV B, the size of the diffusion region scales as l D ∼ (ξ 2 /σ) 1/4 , all other parameters kept constant. With these values the ratio of the diffusion length between colloid-polymer mixtures and molecular liquids is at least two orders of magnitude. Thin films are obtained for high values of both parameters, when advection and viscous stresses are dominant. Our prediction works for both symmetric and asymmetric fluids. From the crossover prediction, we propose that thin films can be assured as long as CaP e > 0.3 for symmetric fluids and CaP e > 0.5 for asymmetric fluids. Anyhow, this values are consistent with the typical experimental regime. For instance, experiments of Ref. [6] were done in a channel of thickness b = 7.95 × 10 −4 m, with a silicone oil with η = 9.3 cP, σ = 20.1 dyn/cm and D = 1.4946 × 10 −6 cm 2 /s [30]. Typical velocities in the experiments ranged from U = 0.01 m/s to U = 0.1 m/s. With these values, we obtain CaP e ≃ 2.6 × 10 4 , which is consistent with our prediction for the thin film regime.
We have examined the steady state of the thin film regime. We have found, in agreement with Ref. [5], that contact line dynamics does not affect the steady state finger shape. The capillary number and the viscosity contrast between the fluids determine the shape of the finger.
For a low-viscosity fluid pushing a high-viscosity one, the finger narrows with increasing capillary number down to a limiting value of 0.57 in units of the channel thickness. These results agree with the numerical results of Ref. [15] for the wide range of capillary numbers considered. Due to computational limitations we do not investigate fingers with very low capillary numbers. Nonetheless, as we recover results from Ref. [15], we expect that the low capillary number limit can be recovered by our method as well.
For fluids with equal viscosities we have found a curve of the finger width as a function of the capillary number that does not follow any previous results. The width of the finger decreases with increasing capillary number, an expected observation, but to a remarkably limiting width of 0.38 in units of the channel thickness. This contrasts with the saturation value of the asymmetric case.
The steady state is independent of whether or not the interface is linearly unstable to the Rayleigh-Taylor instability. This reinforces the conjecture of the steady state being independent of the mechanism that first leads to penetration of one fluid to the other. For low capillary numbers the shape of our fingers is universal and is consistent with the finger shape of Pitts, which suggests that the steady state can be described on simple mechanical equilibrium grounds.
In a future work we will address the problem of viscous fingering allowing for a non-flat leading interface. As we have shown, the shape of the interface across the channel thickness can be controlled by tuning the diffusion strength in the contact line. Hence, it is possible to describe situations in which the leading interface undergoes a fingering process both in the presence and absence of a thin film in the channel thickness direction. In the presence of a thin film, additional control over the interface shape can be gained by choosing between fluids of equal or different viscosities. These features are very convenient to study in detail the three-dimensional effects that arise in the viscous fingering problem.

VI. ACKNOWLEDGMENTS
We acknowledge financial support from Dirección Gen