Three-dimensional aspects of fluid flows in channels. II. Effects of Meniscus and Thin Film regimes on Viscous Fingers

We perform a three-dimensional study of steady state viscous fingers that develop in linear channels. By means of a three-dimensional Lattice-Boltzmann scheme that mimics the full macroscopic equations of motion of the fluid momentum and order parameter, we study the effect of the thickness of the channel in two cases. First, for total displacement of the fluids in the channel thickness direction, we find that the steady state finger is effectively two-dimensional and that previous two-dimensional results can be recovered by taking into account the effect of a curved meniscus across the channel thickness as a contribution to surface stresses. Secondly, when a thin film develops in the channel thickness direction, the finger narrows with increasing channel aspect ratio in agreement with experimental results. The effect of the thin film renders the problem three-dimensional and results deviate from the two-dimensional prediction.


I. INTRODUCTION
Interfacial instabilities in three-dimensional channels give rise to a rich phenomenology in systems that range from nano and microscales[1] to macrometric channels[2, 3,4], and from which a number of practical applications can be drawn.
For instance, controlled drop breakup in microchannels has proved useful in the fabrication of low polydispersity micro-emulsions [5] and in the enhancement of micro-reaction processes [6,7]. In the latter, threedimensional effects are crucial, as they are responsible of a vortex flow structure within the droplet [8] that enhances the mixing process of the reactants.
A widely studied interfacial instability in channels is that of fingering, which occurs whenever a low-viscosity (or high-density) fluid drives a high-viscosity(or lowdensity) one. The instability, first studied by Saffman and Taylor [9], leads to interface dynamics where fingerlike structures emerge and compete. The problem has a steady-state solution, composed by a single finger of constant velocity U and occupies a fraction λ of the width of the channel.
Experimentally, finger growth has been studied mainly in Hele-Shaw cells. These consist of a pair of plates of length L and width W separated by a thickness b. For such systems, it has been pointed out [10] that the stationary finger is determined by a single control parameter, a modified capillary number defined as 1/B = 12Ca/ǫ 2 . For a fluid with viscosity η and surface tension σ, the capillary number, Ca = ηU/σ, measures the competition between driving forces, such as viscous stresses and gravity, and restoring forces, like surface tension. 1/B * Electronic address: rodrigo@ecm.ub.es also includes the degree of asymmetry of the cell, given by the aspect ratio ǫ = b/W . If 1/B is the only control parameter of the system, all experimental data, i.e. all finger widths, should be described by a single curve when plotted as a function of this parameter. Contrary to this view, experiments show that there exists a family of curves λ vs. 1/B for different aspect ratios [11,12]. This fact suggests that a three-dimensional effect, given by the interplay between the dynamics in the channelthickness and in the channel-width, is determinant for the steady-state solution.
Theoretically, fluid-flow in a channel at small velocities pertains to the lubrication regime, in which the flow occurs mainly along the direction of L given that it is much larger than both W and b. Hele-Shaw flows are a limiting case in lubrication theory, where b is much smaller than W . Owing to the smallness of b, the problem is rendered effectively two-dimensional by averaging all fields over the thickness of the channel. Averaging the equations of motion also reduces the interface from a surface to a line, often called the leading interface. In views of the averaged model, three-dimensional effects enter as perturbative corrections to the boundary conditions that hold at the leading interface in terms of Ca and ǫ, particularly to the Gibbs-Thomson condition, which relates the pressure drop across the interface to the interface curvature and surface tension.
Progress towards a three-dimensional description of the problem has been made since the pioneering work of Saffman and Taylor [9], who solved the problem of a stationary finger in the absence of surface tension in two dimensions. McLean and Saffman [13] included the effect of surface tension and were the first to obtain a λ vs. 1/B prediction by solving numerically the twodimensional model. According to their results, λ is a monotonically decreasing function of 1/B that saturates to λ → 1/2 as 1/B → ∞. The prediction of McLean and Saffman is unique in 1/B, so the role of the aspect ratio is precluded from their theory.
The relevance of three-dimensional effects was first suggested by Park and Homsy [14], who pointed out that a thin film of fluid in the channel-thickness direction would contribute to the pressure drop at the leading interface. Using perturbation methods for slightly curved leading interfaces(small ǫ), they found that for low Ca the pressure drop varies as Ca 2/3 , a result that matched the early prediction of Bretherton [15] for capillary tubes.
Sarkar and Jasnow [16] used the modified pressure drop to solve the steady state finger. Their results agreed better with experiments but were restricted to low values of Ca. It was shown by Tabeling and Libchaber [11] that corrections to the pressure drop can be used to reduce three-dimensional experimental data to the twodimensional results of McLean and Saffman. A modified pressure drop can be accounted for as an effective surface tension. Using the correction of Park and Homsy, Tabeling and Libchaber were able to reduce their data to McLean and Saffman results for moderately low values of 1/B, where fitting parameters were used to estimate the correction terms. In an experimental study [12], Tabeling, Zocchi and Libchaber observed that, contrary to the McLean-Saffman prediction, the finger width can go below the one-half limit for sufficiently high 1/B and sufficiently large ǫ.
Reinelt extended the expansion of the pressure drop up to O(1) in Ca and included the effect of larger aspect ratios [17]. Computation of the steady state finger yielded solutions that better agreed with experiments for O(1) values of Ca. For small ǫ Reinelt observed a better agreement between numerics and experiments. However, for relatively large ǫ this agreement is lost.
Higher Ca values have only been explored in the case of flat leading interfaces by Halpern and Gaver [18]. Their numerical results are consistent with results found by Reinelt and Saffman [19] for Ca O(1) and ǫ = 0, and show that the pressure drop is insensitive to the capillary number for Ca > 20.
As an alternative to the sharp interface model, a number of mesoscopic approaches have gained importance in interface dynamics. These are based on order parameter evolution equations of the Cahn-Hilliard type. Being mesoscopic in nature, fluids are separated by diffuse regions instead of sharp interfaces, where the interface boundary conditions arise naturally. All mesoscopic models that address the viscous fingering problem so far are two-dimensional. For fluids of arbitrary viscosities and densities, Folch et al [20,21] used a set of coupled evolution equations for the velocity potential and order parameter that describes accurately the early stages of destabilization of the leading interface, and approaches McLean and Saffman results as the viscosity of the displacing fluid is made negligible. The strict one-sided situation, were one of the fluids is inviscid, was studied by Hernández-Machado et al [22]. They used a single evolution equation for the concentration that includes dynamic effects in the form of chemical potential gradients and described correctly the steady state finger.
In a preceding paper [23], we have shown that a detailed three-dimensional description of fluid-flow in a channel can be done by means of a mesoscopic model which we implement numerically via a Lattice-Boltzmann algorithm. The model considers a fluid-fluid interface in contact with solid boundaries. In contrast to classic approaches, it allows for slip at the contact line by means of a diffusive mechanism inherent to the mesoscopic nature of the interface. This circumvents the complications of contact line dynamics in the classic formulation, associated to the viscous dissipation singularity [24]. In Ref. [23], we focused on the case of a flat leading interface. We showed that depending on the velocity of the contact lines, which we control by modifying the diffusion strength, the interface can either advance as a meniscus or develop as a finger. In the latter case, a thin film of displaced fluid is left adhered to the walls of the channel.
In this paper we extend our Lattice-Boltzmann simulations to the case of a non-flat leading interface, where a viscous finger is expected to appear. Our aim is to provide a detailed description of the mechanisms that affect the steady finger and that cause deviations from two-dimensional results. To do so, we study fingers that form in the meniscus and thin film regimes separately. We cover values of Ca up to O(10) and explore various aspect ratios.
The paper is organized in the following manner. In Sec. II we present the governing equations of the system which we solve numerically by means of a Lattice Boltzmann algorithm, presented in the preceding paper [23]. Results are presented in Sec. III. In Sec. III A we describe the simulation strategy and parameter steering procedure. As a validation test, in Sec. III B we compute the dispersion relation of the interface in the twodimensional limit and compare it to the analytic prediction of the Saffman-Taylor problem. Sec. III C is devoted to the study of stationary viscous fingers; in Sec III C 1 we focus on fingers pertaining to the meniscus regime in the channel thickness, which we found to be effectively twodimensional, while in Sec. III C 2 fingers in the thin film regime are studied. We find that fingers in the thin-film regime are three-dimensional and cannot be described by the two-dimensional theory in general. A discussion of our results where we compare with previous experiments is presented in Sec.IV. In Sec. V we present the conclusions of this work.

II. GOVERNING EQUATIONS
We consider the motion of two viscous fluids, whose dynamics are governed by the Navier-Stokes equations, Here v is the fluid velocity, P is the pressure, ρ is the density, η is the fluid viscosity and g is the acceleration of gravity. The extra term, φ ∇µ, is mesoscopic and accounts for interfacial forces. Here, φ( r, t) is an order parameter and µ(φ) is the chemical potential. φ has the property of being uniform in the volume of each phase and non-uniform in an interfacial region of typical size ξ. In the present case, volume values are chosen as φ = ±1 for the displacing and the displaced fluid respectively, so the interface is located at φ = 0. The size of the interface is set to ξ = 0.57.
The dynamics of φ obey a convection-diffusion equation, where M is a mobility coefficient. In equilibrium, the pressure and chemical potential minimize a free energy functional, from which explicit expressions P (ρ) and µ(φ) can be derived. For further details the reader is referred to Ref. [23]. We work on a linear channel, composed by two solid plates of width W and length L parallel to the xy plane, separated by a distance b, as depicted in Fig. 1. There exist two principal directions in the system: a lateral direction, parallel to the xy plane, and a transverse one parallel to the xz plane. We will denote these by subscripts and ⊥ respectively.
The impenetrability and stick boundary conditions at the walls are enforced by setting v(x, y, As for the fluid-fluid boundary, the Gibbs-Thomson relation is recovered by integrating Eq.(1) across the interfacial region, where σ is the surface tension and R α is the radius of curvature of the interface in the direction α.
We now briefly review the classic treatment of the problem. First, v ⊥ is assumed to be much smaller than v , which in turn is expected to be parabolic in z. As a result, Eq.(1) is recast in the form of an average velocity field which holds in the volume of each fluid, called where triangular brackets denote an average over the channel thickness. Under these conditions, R ⊥ is expected to be much larger than R . Hence, in the twodimensional theory the Gibbs-Thomson relation is simplified to ∆P = σ/R . Corrections to this expression arise whenever 1/R ⊥ is not negligible.
For such cases, Libchaber and Tabeling [11] have proposed that thin film effects can be accounted for by defining an effective surface tension so the two-dimensional form of the Gibbs-Thomson condition is recovered. For this purpose, they used the estimation of Park and Homsy [14] of the pressure drop for Ca → 0 and slightly curved leading interfaces(ǫ → 0), As a result, their experimental results collapsed to the McLean-Saffman curve when using the corresponding definition of the control parameter, 1/B * = (σ/σ * )1/B. We solve numerically Eqs.
(1) and (2) by means of a Lattice-Boltzmann algorithm. For further details of the method, the reader is referred to the preceding paper [23].

A. Simulation Parameters and Setup
The traditional description of the viscous fingering problem corresponds to situations in which the relevant forces at play are viscous stresses and capillarity. For the particular case of fingering in a Hele-Shaw cell these forces are expressed in terms of a modified capillary number [10] 1/B = 12(W/b) 2 (∆ηU +∆ρgb 2 /12)/σ, where ∆η and ∆ρ are the differences in viscosity and density between the fluids.
To ensure that capillarity and viscous forces dominate the dynamics of the fluids, inertia must be small compared to both of these forces. We enforce this situation by neglecting the convective term in Eq. (1). As for compressibility, we consider low Mach number flows, which we achieve by keeping U ≪ c s . For our scheme, it suffices to set U ≤ 0.01.
Our goal is to explore the viscous fingering problem for a wide range in 1/B. Due to computation resource limitations, ǫ is restricted to O(10) at most for the majority of runs. To achieve high values of 1/B, say O(10 3 ), Ca The channel is implemented as follows. We set a rectangular simulation box of dimensions N x × N y × N z. Due to the flow symmetry, we simulate only one fourth of the real channel by setting boundary conditions as follows: ∂ y ρv y (x, y = 0, z) = ∂ y ρv y (x, y = W/2, z) = 0, ∂ y φv y (x, y = 0, z) = ∂ y φv y (x, y = W/2, z)0, ∂ z ρv z (x, y, z = 0) = ∂ z φv z (x, y, z = 0) = 0.

B. Linear Stability in the two-dimensional limit
We first verify the linear stability of the interface, i.e., the behavior of an initially flat interface that has been subjected to a small perturbation. We study fluids of equal viscosities, so the instability is gravitationally driven. This is done by fixing the body force term in Eq. (1) as ρg = 8η/b 2 U exp (φ + 1)/2, where U exp is the maximum expected velocity for a Poiseuille flow. The modified capillary number reduces to 1/B = W 2 ∆ρg/σ, In this case, the linear stability analysis of the interface evolution of the averaged equations yields the dispersion relation where ω is the exponential growth rate of a sinusoidal perturbation to the flat interface solution. The perturbation is characterized by its wavelength, Λ = 2π/k. By considering dimensionless frequencies ω ′ = ω/(U/2W )B 1/2 and wavenumbers, k ′ = W B 1/2 k, the dispersion relation becomes universal, i.e., ω ′ (k ′ ) = k ′ (1 − k ′ 2 ). We prepare a base flow corresponding to a flat interface in the xy plane that propagates at constant velocity. The interface in the xz plane is nearly flat throughout the simulation, so the system is effectively two-dimensional. Once the base flow is fully developed, the interface is shifted according to a single mode perturbation of wavelength Λ = W and an initial small amplitude. We follow the evolution of the amplitude, A(t), which is measured as A(t) = x tip (t)−x(t), where x tip andx(t) are the interface tip and mean interface positions respectively. The growth rate, ω, is extracted as a linear fit of log(A(t)) vs t. Fig. 2 shows a comparison between the universal dispersion relation and our results. To quantify the degree of accuracy of these results, we fit our data to the general form ak ′ (b − ck ′ 2 ). We find a most unstable mode at k ′ max ≃ 0.56 and a first unstable mode at k ′ 0 ≃ 0.96, both 4% below the exact result.

C. Viscous Fingers
In a preceding study [23], we have shown that it is possible to control the generation of a thin film in the channel by adjusting the diffusivity of the order parameter. Although for usual experimental conditions this is not a relevant parameter (it might be relevant in nano-channels), it gives the possibility of elucidating the role of a thin film in viscous fingers. Diffusivity is accounted for by a Péclet number, P e = U b/D, where D is the diffusion coefficient. By combining the effects of P e and Ca, one can either suppress or induce the formation of a thin film. In particular, a small value of the product CaP e results in the suppression of thin films, while the contrary is obtained for high CaP e. Results from the preceding work give a penetration threshold of CaP e ≃ 10 −1 .
The strategy is to study first fingers for which CaP e ≤ 10 −1 and then extend this simulations to CaP e ≫ 10 −1 .

Meniscus Regime
We first study fingers for which no film of displaced fluid develops in the xz plane of the channel. We carry out simulations with modified capillary numbers in the   Table I. The error (small bar at the right bottom) corresponds to one lattice spacing. The larger bar indicates the size of the diffuse interface, approximately 3ξ.
range 100 ≤ 1/B ≤ 6000. We have studied different geometries, ranging from ǫ = 0.17 to ǫ = 0.04. The aspect ratio is decreased by decreasing the channel thickness. We summarize the simulation parameters in Table I. For each run we observe the usual phenomenology for the leading interface. During the early stages of interface evolution, the amplitude of the perturbation grows until a finger emerges and widens. This stage is followed by a relaxation of the interface shape, until a Saffman-Taylor finger develops. The finger propagates with a steady velocity U , leaving behind a growing region where the finger has flat sides. In this region a constant finger width, λW , can be defined. As for the channel thickness, we observe that the initially flat interface rapidly relaxes to a meniscus, which also has a steady shape. In Fig. (3) we show a three-dimensional plot of the interface for run (a) in Table I at two different times. In the plot we show both the contact lines and the leading interface; both contact lines follow the leading interface.
To check for consistency in the steady state solution we use the semiempirical interface profile obtained by Pitts [25], which reproduces experimental results accurately for a wide range of finger widths. The equation for the interface shape reads, where x ′ and y ′ measure the distance from the finger tip in units of half the channel width. The natural scalings in this equation are πx ′ /2λ and πy ′ /2λ. Consequently, all profiles should collapse into a single curve if these scalings are used. Fig. 4 shows plot of interface profiles corresponding to runs of Table I. As expected, all interface profiles fall in the same universal curve within error. In addition, our collapse is in fair agreement with Eq. (7). The selection rule in the viscous fingering problem is expressed as the functional dependence of the finger width with the modified control parameter. We compare our results with the λ vs. 1/B prediction of McLean and Saffman. We find that runs with ǫ = 0.17 show wider fingers than predicted, while runs with smaller ǫ agree better with the two-dimensional result. Even in the absence of a thin film, the xz interface projection has a certain curvature. This can be accounted for by defining an effective surface tension in terms of the radii of curvature of the interface, which we are able to measure directly. The effective surface tension then reads σ * = σ(1 + R /R ⊥ ). The correction factor in this expression is given by the quantity in parentheses, which increases for strongly curved meniscus. The rescaled control parameter then reads 1/B * = (1/B)/(1 + R /R ⊥ ). Of course this correction should be more evident in the low 1/B region, where λ varies rapidly with the modified control parameter. In Fig. 5 we show a plot of λ vs. 1/B * . Points fall on the McLean-Saffman curve for the wide range of 1/B * considered, regardless of the aspect ratio.

Thin Film Regime
We now extend our simulations to fingers in the film regime. Penetration in the xz plane occurs for high CaP e, so we choose to sample 1/B at fixed D. Conse- quently, CaP e increases with increasing 1/B. To resolve the thin film correctly we must take into account the finite size of the interface. As explained in Ref. [23], the thin film is insensitive to the channel thickness already for b = 23. We therefore choose sufficiently thick channels. We explore a wide range of aspect ratios, 0.25 ≤ ǫ ≤ 1.0 and modified capillary numbers, 800 ≤ 1/B ≤ 5300. We first explore the CaP e ∼ O(1) region, close to the penetration threshold. In Fig. 6 we show interface projections in the xy and xz planes located at z = b/2 and y = W/2 respectively. We show two sets of interfaces, corresponding to two different CaP e values; (a)CaP e = 0.85 and (b)CaP e = 4.44. In Fig. 6(a) the interface in the xz plane presents a penetrating structure, but a well developed film is absent. The finger in the xy plane is not well developed either, and it presents an anomalous tip. Conversely, in Fig. 6(b) both interface projections describe well developed fingers. It is then clear that deviations from the Saffman-Taylor finger in the xy plane are correlated to the interface structure in the xz plane. An interesting feature of the run corresponding to Fig. 6(a) is that the the xz interface structure is persistent. This means that the length of the finger in the xz plane is constant in time, a consequence of the slip velocity of the contact line. The diffusion strength is not large enough to maintain a meniscus, which in the one hand makes the slip velocity smaller than the channel velocity. Nevertheless, as the interface relaxes to a thin film shape, curvature deviations from equilibrium increase the slip velocity, making the contact line advance to restore the meniscus shape.
We next explore the range CaP e ≥ O(10) for which simulation parameters and observed finger widths are summarized in Table II. In Fig. 7 we present snapshots of the three dimensional interface at two different times for run (b) in Table II. The first snapshot corresponds to the early stage of the finger formation. Looking at the interface projections in the xy plane, we see that the contact line(light line) is close to the leading interface(dark line) and no film is present in the xz plane. In the next snapshot the contact line has moved away from the tip, thus giving rise to the growth of a wetting film. The shape of the finger is in agreement with the typical morphology found in experiments. To illustrate this, in Fig. 8 we compare the shape of the finger to Eq. (7). Within error, our profiles are consistent with Pitts shape. Fig. 9 shows the measured finger width as a function of 1/B. The lowest aspect ratio we have considered corresponds to ǫ = 0.25(runs (a)-(d) in Table II). We see that for all 1/B values considered the finger width falls above the McLean-Saffman prediction. We increase the aspect ratio to ǫ = 0.35(run (e) in Table II). As a result, the measured finger width decreases. Runs for which ǫ is larger confirm this tendency in a systematic way. Tests (f)-(j) in the same table correspond to a fixed value of 1/B with increasing ǫ. We find that for sufficiently large ǫ the finger width goes below the one-half theoretical limit of McLean and Saffman.  Table II. The bars in the bottom at the right indicate the error bar and diffuse interface size as in Fig. 4.

IV. DISCUSSION
Our results show that the finger width decreases with increasing aspect ratio. To maintain 1/B fixed while  varying the aspect ratio of the channel, one has to vary Ca accordingly. As a consequence, the film thickness and the capillary pressure are altered. If we increase the aspect ratio(as in the high-1/B region in Fig. 9), then Ca must decrease to keep 1/B fixed. As a consequence, the film thickness and the capillary pressure decrease as ǫ increases, which is consistent with a narrower finger. This behavior has been observed, for instance, in experiments by Tabeling, Zocchi and Libchaber [12], and addressed in numerical calculations by Reinelt[17] where the effect of the thin film was introduced perturbatively in the two-dimensional model. Experiments suggest that, for high 1/B, increasing the cell aspect ratio has a thinning effect on the finger, which is what we observe in our simulations. Results of Reinelt suggest the opposite tendency.
The aforementioned experiments were done at small ǫ and Ca, and at high viscosity contrast, defined as c = (η 2 − η 1 )/(η 2 + η 1 ). As we have shown in Ref. [23], the thin film thickens as c → 1. Under these conditions, experiments show that the finger width goes below the one-half limit even for cells with ǫ = 0.009. This is due to the small thickness of the film which is a consequence of the low Ca and high c values used in experiments. In our simulations the thin film is about t/b ≃ 0.25, much thicker than the experimental value, t/b ≃ 0.05. As a consequence, curvature effects in our simulations are strong enough to keep the finger width above one half even for high values of ǫ. To achieve the experimental regime thinner film should be considered. We have considered a cell aspect ratio of ǫ ≃ 0.05 and c = 0.9. Nevertheless, Ca is still too large, the film is then thick enough to keep us in the low 1/B * regime, where the finger width is still larger than one half of the channel width. Due to computational limitations we do not explore smaller ǫ.
The fact that for a given 1/B there exist different finger widths for different aspect ratios raises the doubt of 1/B as being the only control parameter present in the system. To this end, we compute the rescaled surface tension σ * = σ 1 + R /R ⊥ , where the radii of curvature are measured at the finger tip. We then rescale the control parameter according to 1/B * = (σ/σ * )1/B. In Fig. 10 we show a plot of the finger width as a function of the rescaled control parameter. At low 1/B * , results agree with McLean-Saffman results. We conclude that in this region the finger is effectively two dimensional and that three-dimensional effects can be accounted for even at Ca ∼ 1.
At high values of the rescaled control parameter, our results deviate systematically from the McLean-Saffman curve, until the finger width goes below the one half limit of the two-dimensional theory. This behavior is qualitatively different from the one found for the meniscus regime, in which the McLean-Saffman curve could be recovered at any value of 1/B * . Hence, we conclude that deviations from two-dimensionality are caused by the thin film.
An important feature in the λ vs. 1/B * plot is that finger width appears to be determined by 1/B * uniquely. This suggests that 1/B * is the only control parameter of the problem.
We have explored a region of values of the aspect ratio between the Saffman-Taylor(ǫ → 0) and Rayleigh-Taylor(ǫ = 1) limits of the fingering instability. In both limits, the relevant control parameter appears to be an effective modified capillary number. In addition, the interface shape is remarkably universal, as suggested by Figs.4 and 8.

V. CONCLUSIONS
We have carried out a detailed study of the viscous fingering problem in three-dimensional channels for fluids of different densities and viscosities. We have studied the single finger solution for systems in which either a thin film develops across the channel thickness or a meniscus is stabilized.
For systems in which no thin film is present, McLean-Saffman two-dimensional results describe the dependency of the finger width as a function of a rescaled modified capillary number, 1/B * , which has a correction that depends curvature of the interface direction of the channel thickness. This holds for arbitrary high values of 1/B * , evidencing that a complete displacement across the channel thickness renders the problem two-dimensional.
We have extended our studies to situations where a thin film develops across the channel. We find different values of the finger width when changing the channel aspect ratios at fixed modified capillary number, an observation that is consistent with previous experiments [12]. This non-uniqueness seems to disappear as the control parameter is corrected by curvature effects associated to the thin film, i.e., when the finger width is compared to 1/B * .
For low 1/B * , the finger width collapses to the McLean-Saffman curve. However, at high 1/B * the finger width deviates from this curve, and goes bellow the limit of one half in units of the channel width.
Our work is done at high values of the capillary number. Consequently, the effective capillary pressure in our simulations is large enough to keep the finger above the one-half limit of the two-dimensional theory for high values of the channel aspect ratio. Experiments in Refs. [11,12] were done in cells with ǫ ≃ 0.03 and at Ca ≃ 10 −3 , a regime that falls beyond the scope of this work for computational reasons. Nonetheless, for low 1/B * , we recover the same results, indicating that the same mechanisms hold, even if the actual aspect ratio and capillary number are not the same in experiments and simulations.
To our knowledge, experiments of fingering in high aspect ratio channels have not been conducted so far. Our results could be confirmed, for instance, in microchannels, where the aspect ratio is typically large and in which the meniscus to thin film transition could be observed. This is then an open question for experimentalists to confirm.