Pore-filling events in single junction micro-models with corresponding lattice Boltzmann simulations

The aim of this work is to better understand fluid displacement mechanisms at the pore scale in relation to capillary-filling rules. Using specifically designed micro-models we investigate the role of pore body shape on fluid displacement during drainage and imbibition via quasi-static and spontaneous experiments at ambient conditions. The experimental results are directly compared to lattice Boltzmann (LB) simulations. The critical pore-filling pressures for the quasi-static experiments agree well with those predicted by the Young–Laplace equation and follow the expected filling events. However, the spontaneous imbibition experimental results differ from those predicted by the Young–Laplace equation; instead of entering the narrowest available downstream throat the wetting phase enters an adjacent throat first. Thus, pore geometry plays a vital role as it becomes the main deciding factor in the displacement pathways. Current pore network models used to predict displacement at the field scale may need to be revised as they currently use the filling rules proposed by Lenormand et al. (J. Fluid Mech., vol. 135, 1983, pp. 337–353). Energy balance arguments are particularly insightful in understanding the aspects affecting capillary-filling rules. Moreover, simulation results on spontaneous imbibition, in excellent agreement with theoretical predictions, reveal that the capillary number itself is not sufficient to characterise the two phase flow. The Ohnesorge number, which gives the relative importance of viscous forces over inertial and capillary forces, is required to fully describe the fluid flow, along with the viscosity ratio.

FIGURE 1. Schematic drainage/imbibition capillary pressure curve. During primary drainage the capillary pressure P c systematically increases with increasing non-wetting phase saturation to connate water saturation (S WC ), whereas the capillary pressure decreases with increasing wetting phase saturation during imbibition. Lenormand et al. (1983) illustrated that the Young-Laplace equation was sufficient to describe all menisci displacement mechanisms. The Young-Laplace equation (Rowlinson & Widom 1982) states where P c is the capillary pressure, P nw and P w are the non-wetting phase (NWP) and wetting phase (WP) pressures respectively, γ is the surface tension and r 1 , r 2 are the principal radii of curvature of the fluid interface. For throats with a rectangular cross-section this becomes P c = 2γ cos θ(1/d + 1/w t ), where θ is the contact angle and d, w t are the depth and width of the throat. For primary drainage the porous medium is initially saturated with WP before NWP is forcibly injected, displacing the WP and thus causing the capillary pressure to systematically increase with increasing NWP saturation, see figure 1. Considering a pore junction with different throat widths, where the NWP enters from one of the throats, the Young-Laplace law (1.1) says that the NWP should select the downstream throat with the lowest capillary entry pressure first, i.e. by entering the widest throat. In the case of CO 2 sequestration in a saline aquifer, once injection of the NWP (CO 2 ) (Chiquet & Broseta 2007;Doughty, Freifeld & Trautz 2008;Hesse, Orr & Tchelepi 2008;Chalbaud et al. 2009) has stopped, the WP (brine) re-enters the porous medium via imbibition. During the imbibition process, the capillary pressure systematically decreases with increasing WP saturation, see figure 1. This means that, following (1.1), we expect that the WP will enter the narrowest downstream throat first, before sequentially filling all other throats in order of increasing diameter. In reality, displacement processes are more complex as they are dependent on the pore geometry along with the WP and NWP location. In order to observe these displacement mechanisms experimentally, we used: (i) micro-fluidic devices as their transparent nature makes visualisation of fluid dynamics relatively simple and (ii) single pore geometries, rather than more complex networks, as single pore geometries allow us to accurately control the displacement processes in each throat individually; see figure 2 for schematic micro-model designs used. Furthermore, the models' characteristics, like pore geometry, surface energy and roughness, can be carefully controlled. The above enabled us to examine real-time primary drainage and imbibition events systematically via both quasi-static and spontaneous methods. The quasi-static experiments are performed in a sequence of pressure steps, with hydrostatic equilibrium reached before progressing to the next pressure step. These quasi-static experiments are relevant to displacement processes occurring in the far field during CO 2 storage and enhanced oil recovery operations, associated with low Reynolds numbers. Displacement processes within pore network models (Øren, Bakke & Arntzen 1998;Øren & Bakke 2003;Valvatne & Blunt 2003;Sorbie & Skauge 2012) are performed in this fashion. The dynamic experiments were conducted as they are relevant to drainage and imbibition processes occurring near the well bore at relatively high Reynolds numbers. Further details on the experimental methods will be provided in the next section.
To further elucidate our experimental findings we compare experiments to lattice Boltzmann simulations in the same geometries. These provide the opportunity for a detailed investigation of the problem by varying the parameters under investigation over a wider range, not easily accessible to the experiments. The overall goal of this study is to challenge the pore-filling rules on which network models are based by comparing the sequence of throat filling in simple geometries for which the Laplace pressures can be calculated exactly and without ambiguity. This paper is organised as follows; in the next section we provide the details of the experiments and the numerical scheme. We present and discuss our results in § 3. Finally conclusions drawn from this work are discussed in § 4. To conduct these experiments we had specifically designed micro-models fabricated in poly(methyl methacrylate) (PMMA), figure 2(a-c), by Epigem (Redcar). The fabrication procedure involves defining the pattern in the base layer of the model by using SU-8 (an epoxy based photoresist) via photolithography, which is achieved in two stages. Initially an under-layer is deposited (spin coated then dried) and fully cured before a secondary layer is deposited in the same way. The pattern is then created by exposing the coated models through a photo-mask and developing it to form the features. The top layer has a partially cured SU-8 layer deposited on the underside (this will form the top of the pattern) before the inlet/outlet holes are drilled into it. Finally, the base and top layer are assembled; figure 2(e) illustrates the cross-section of an assembled model. All the models have been chemically treated by Epigem, in order to increase the hydrophobicity of the surface. Additionally all models have an approximate etch depth of d = 50 µm. The designs are intended to explore pore geometry and the influence of varying throat diameters, providing a range of different capillary entry pressures.

Experimental set-up
All of the fluid displacement observations were captured via a high-speed video microscope (FastCam MC2.1, Photron) which is housed within a laminar flow cabinet (PurAir-48, Air Science). This is necessary due to the intricate nature of the micro-models, where unnecessary exposure to dust can lead to blockages rendering the micro-models unusable. Due to the hydrophobic nature of the micro-models, n-decane (viscosity η w = 0.85 mPa s (Dymond & Øye 1994), surface tension γ = 0.024 N m −1 (Kuhn, Försterling & Waldeck 2009), density ρ = 730 kg m −3 , 99 %, Sigma-Aldrich) and air (viscosity η nw = 18.2 µPa s, density ρ = 1.2 kg m −3 ) were selected as the WP and NWP respectively. However, due to the unique fabrication of each micro-model, the etch depth/width and contact angle varied for each chip, with the contact angle always measured through the denser phase (Lyons 2009), resulting in different capillary entry pressures. All experiments were conducted at ambient conditions.

Experimental procedure -primary drainage and imbibition
The primary drainage experiments began with the models fully saturated with the WP. The NWP entered the model via either: (i) quasi-static displacement, achieved by gradually decreasing the WP pressure by siphoning the WP into a reservoir located below the model via the narrowest throat, or (ii) dynamic displacement -injection of the NWP at a constant flow rate of 0.5 µl min −1 via a programmable syringe pump (BS-8000, Braintree Scientific Ltd).
For imbibition the models were initially saturated with the NWP. During quasi-static displacement the model was connected to a reservoir of the WP and P w was then gradually increased by raising the height of the reservoir. In contrast, spontaneous displacement was attained by placing a droplet of the WP over an inlet port which then spontaneously penetrated into the model.

Numerical method
In this section we describe the numerical method we shall use, starting with the thermodynamics in § 2.2.1, the dynamical equations of motion in § 2.2.2 and the lattice Boltzmann implementation in § 2.2.3.

Thermodynamics of the fluid
The equilibrium properties of a binary fluid can be described by a Landau free energy functional (Briant & Yeomans 2004) (2.1) The first term in the integrand is the bulk free energy density given by where φ is the concentration or order parameter, ρ is the fluid mass density and c is a lattice velocity parameter. A is a constant with dimensions of energy per unit volume.
This choice of f b allows binary phase separation into two phases with bulk equilibrium solutions φ eq = ±1. The position of the interface is chosen to be the locus φ = 0. The term in ρ controls the compressibility of the fluid (Kendon et al. 2001).
The presence of interfaces is accounted for by the gradient term (κ φ /2)(∂ α φ) 2 , which penalises spatial variations of the order parameter φ. This gives rise to the interface tension γ = 8κ φ A/9 and to the interface width ξ = κ φ /A (Briant & Yeomans 2004).
The final term in the free energy functional, equation (2.1), describes the interactions between the fluid and the solid surface. Following Cahn (1977), the surface energy density is taken to be of the form f s = −hφ s , where φ s is the value of the order parameter at the surface. Minimisation of the free energy gives an equilibrium wetting boundary condition (Briant & Yeomans 2004) 3) The value of the parameter h (the surface excess chemical potential) is related to the equilibrium contact angle θ eq via ( where α = arccos(sin 2 θ eq ) and the function sign returns the sign of its argument. This choice of the free energy leads to the chemical potential (2.5) and the pressure tensor (Anderson, McFadden & Wheeler 1998) where p b = (c 2 /3)ρ − (Aφ 2 )/2 + (3/4)Aφ 4 is the bulk pressure.

Equations of motion
The hydrodynamic equations for the system are the continuity, (2.7), and the Navier-Stokes, (2.8), equations for a non-ideal fluid (2.8) where u, P, η are the fluid velocity, pressure tensor and dynamic viscosity respectively. For a binary fluid the equations of motion are coupled with a convection-diffusion equation, that describes the dynamics of the order parameter φ. M is a mobility coefficient.

Lattice Boltzmann method
The equations of motion are solved using a standard free energy lattice Boltzmann algorithm for a binary fluid (Briant & Yeomans 2004). In particular we use a three-dimensional model with 19 discrete velocity vectors (D3Q19) and adopt a multiple relaxation time (MRT) (D'Humières et al. 2002) approach for the evolution of the distribution functions, f i , associated with the fluid density ρ. Following Pooley, Kusumaatmaja & Yeomans (2009), the relaxation times responsible for generating the viscous terms in the Navier-Stokes equation are set to τ f (based on the fluid viscosity η = ρ(τ f − 0.5)/3), those related to conserved quantities to infinity and all the others, which correspond to non-hydrodynamic modes, to unity. A single relaxation time approach is sufficient for the order parameter φ. The relaxation time for the evolution of the distribution functions, g i , associated with the concentration φ is set to τ g = 1. As shown by Pooley et al. (2009), this approach suppresses spurious currents at the contact line, while improving significantly the numerical stability of the algorithm as well (Lallemand & Luo 2000). For a detailed description of the lattice Boltzmann (LB) implementation we refer the reader to (Briant & Yeomans 2004;Yeomans 2006;Pooley et al. 2008Pooley et al. , 2009).
Numerically we solve the equations of motion using graphics processing units (NVIDIA Tesla K40 GPU cards), taking advantage of the fact that the LB method is particularly well suited for computations on a parallel architecture (Gray, Cen & Boek 2016). Running the simulations in parallel on 8 Tesla K40 cards reduces the required computational time to a few days for the most computationally intensive simulation.

Results and discussion
In this section we will present our results. We start with primary drainage, which also serves as a validation for our experiments. Then we turn our attention to the imbibition case, where we examine pore filling events and investigate the role of pore geometry in determining the displacement pathways.

Primary drainage
The model was initially filled with the WP. For the quasi-static displacement, the P w was gradually decreased, allowing the NWP to enter the largest throat first, before sequentially displacing the WP from the throats in decreasing size order, as can be seen in figure 3. This type of displacement is referred to as 'piston-type' motion Pore-filling events in single junction micro-models 557 FIGURE 3. Quasi-static primary drainage images in a micro-model (Geometry 2) with etch depth d = 45 µm and throat widths w t : 33, 51, 65, 107 µm (WP-white, NWP-grey). Contact angle 26 • .

Throat
Young-Laplace equation Young-Laplace Experimental pressure (Pa) Critical pressure (Pa)  (Lenormand et al. 1983), when the NWP enters the throat filled with the WP only if the capillary pressure is equal to or greater than a given value P p = P nw − P w = 2γ cos θ eq (1/d + 1/w t ), which we call 'the critical pressure'. The calculated theoretical and experimental critical pressures for drainage within the throats are displayed in table 1. Generally there is good agreement between the two. For the dynamic drainage, the NWP is injected into the chip at a constant flow rate of 0.5 µl min −1 . This leads to the displacement of the WP at a constant mean velocity u in each throat, prior to and after the square pore body. The corresponding capillary number for the interface motion in the throat prior to the pore body is Ca = η w u/γ ∼ 10 −4 . Typical values for the ratio of viscous to capillary forces at the pore scale, quantified by the capillary number, are in the range of 10 −3 -10 −10 , depending on the distance from the injection point in the well bore (Blunt & Scher 1995). Again, as can be seen in figure 4(a), the NWP displaces the WP via the largest downstream throat first, as predicted by (1.1). To validate the displacement sequence, we have carried out LB simulations in exactly the same geometry, the same value for Ca and contact angle as in the experiment. The fluid flow was driven by applying velocity boundary conditions (Hecht & Harting 2010) at the inlet and outlet of the simulation domain to match the experimental conditions. The results are presented in figure 4(b) and display a good agreement with the experimental displacement sequence.

Quasi-static imbibition
Pore-filling events are dependent on the number of throats connected to the pore that are occupied with the NWP along with their orientation. Generally, pore-filling events are designated as I 1 , I 2 , I 3 . . . I n , with n indicating the number of throats 558 I. Zacharoudiou, E. M. Chapman, E. S. Boek and J. P. Crawshaw Pore-filling Young-Laplace Young-Laplace Experimental pressure (Pa) events for pore-filling event Critical pressure (Pa) TABLE 2. Critical pressures for each pore-filling event in a micro-model (Geometry 1) with equal throat widths. Contact angle: 28 • , etch depth d = 56 µm, throat width w t = 73 µm, pore width: 238 µm, interfacial tension: 0.024 N m −1 . The critical radii in the plane of the micro-model r 2 = r used in the calculation of the critical pressures is shown in figure 5. The principal radii of curvature, equation (1.1), are r 1 = d/2 cos θ eq , r 2 = r. Experimental error of ±7 Pa is attained by the 1 mm accuracy the height of the reservoir can be set to.
filled with the NWP (Lenormand et al. 1983). Here I 1,2,2A,3 pore-filling events were investigated via quasi-static imbibition experiments, by gradually increasing P w , in a micro-model with equal throat widths (Geometry 1: etch depth d = 56 µm, throat width w t = 73 µm, pore width 238 µm), see figure 5. At a critical capillary pressure the WP spontaneously displaces the NWP via 'piston-type' displacement. In all cases the critical pressure of the pore-filling event calculated from the Young-Laplace equation is very close to the experimentally determined critical pressure, as shown in table 2. Good agreement is achieved as we know the exact pore geometry and can observe the meniscus shape, allowing the critical radii to be well defined, as illustrated for each case in figure 5. Quasi-static experiments were also conducted in a micro-model (Geometry 3) with unequal throats (etch depth d = 55 µm, throat widths w t = 27, 47, 60, 100 µm pore: 155, 236 µm). The predicted displacement sequence for this geometry, . (Colour online) Quasi-static I 1 , I 2 , I 2A and I 3 imbibition experimental results (WP -white, NWP -grey). The critical radii in the plane of the micro-model, r 2 = r, used to calculate the displacement pressures for each case are also illustrated together with r 1 = d/2 cos θ eq . using the Young-Laplace equation and a contact angle of 16 • , is first snap off (P s ) in the narrowest throat at 1240 Pa, followed by pore body filling via 'piston-type' displacement (P p ) at 869 Pa, see figure 6(a). The snap-off pressure, P s = 2γ (cos θ eq − sin θ eq )/ min(d, w t ), was estimated from the pressure at which two growing corner menisci meet on the channel wall (Valvatne & Blunt 2004). Snap off is able to occur during the quasi-static experiments, as there is time for the advancing WP films to develop ahead of the main meniscus, which swell and become unstable (Bico & Quéré 2003;Kavehpour et al. 2003), see figure 6(b), and because 'piston-type' displacement is not possible for topological reasons (Lenormand et al. 1983). The WP films can be clearly seen in figure 7 as the darker lines outlining the model geometry.

Spontaneous imbibition
Extending our research to spontaneous imbibition, we investigated the I 3 displacement mechanism, which was compared to LB simulations, figure 8. During our spontaneous imbibition experiments we observed that the advancing meniscus did travel down an adjacent throat (T 2 ) first instead of filling the smallest throat, as the WP films do not have time to develop ahead of the advancing meniscus. Indeed, Lenormand (1990) noted that when no wetting films are present along the corners of the models, the imbibing WP should enter an adjacent throat first, but no direct evidence was provided. Here we confirm this hypothesis directly and show our results in a series of snapshots from experiment and corresponding LB simulations in figure 8.    Thus, displacement does not occur in the order predicted by the Young-Laplace equation. This is understandable as all the downstream throats have lower critical entry pressures than the pore body. Therefore, as soon as the critical pore-filling pressure has been exceeded, the WP can enter any of the downstream throats. Hence, in the case of spontaneous imbibition, the pore body geometry becomes a key factor in determining in which order the downstream throats will fill.
Capillary-filling dynamics in a rectangular throat A reasonable assumption then could be that the dynamics of imbibition in the throat prior to the junction (T 3 ) may affect the pore filling sequence and the selection of the displacement pathway. Hence, we turn our attention to the imbibition dynamics. We recently (Zacharoudiou & Boek 2016) examined the dynamics of capillary filling in two-dimensional channels and covered both: (i) the limit of long times for both high and low viscosity ratios r η = η w /η nw and (ii) the limit of short times, demonstrating that the free energy LB method can capture the correct dynamics for the process. We recall that in the limit of high viscosity ratios and long times, when the total time is much larger than the viscous time scale t v ∼ ρL 2 s /η w (Quéré 1997;Stange, Dreyer & Rath 2003), the Lucas-Washburn regime (l ∼ sqrt(t)) (Lucas 1918;Washburn 1921) is expected. In the limit of short times two regimes, namely (i) inertial regime (l ∼ t 2 ) and (ii) visco-inertial regime (l ∼ t) can precede the Lucas-Washburn regime (Dreyer et al. 1994;Stange et al. 2003). The hydraulic diameter D h = 2dw t /(d + w t ) is used as the characteristic length scale L s for evaluating t v . Ichikawa, Hosokawa & Maeda (2004) examined the interface motion driven by capillary action in the case of three-dimensional channels with a rectangular cross-section. Using the analytical solution of Brody, Yager & Austin (1996) for the velocity profile in a rectangular channel of aspect ratio = d/w t and balancing the relevant forces, they estimated the dimensionless relation for the interface position as a function of time l * 2 = cos θ eq (t * − 1 + e −t * ). (3.1) The rescaled time and length are defined as t * = t/t c and l * = l/l c respectively, where Here, we first compare our results with the theoretical prediction of Ichikawa et al. (2004), equation (3.1), before examining whether the capillary-filling dynamics affects the displacement pathway after the junction. We considered equilibrium contact angles of θ eq = 30 • and 16 • (experimental condition). Starting with a situation with a high viscosity ratio r η = η w /η nw = 500 by choosing relaxation rates τ f ,w = 1.5 and τ f ,nw = 0.502, we decrease the viscosity ratio to r η = 5 by decreasing η w , while η nw is kept fixed. This choice decreases the rate of viscous dissipation in the WP as the viscosity ratio decreases from r η = 500 to r η = 5, allowing for more energy to be available to the interface motion as it enters the junction region. Figure 9(a) shows the length of the penetrating WP column in throat T 3 for simulations with varying viscosity ratio r η and θ eq = 30 • . In the limit of long times and high viscosity ratios the Lucas-Washburn regime (l ∼ sqrt(t)) (Lucas 1918;Washburn 1921) is clearly observed. The viscous time scale t v is of the order of 10 4 in lattice units (l.u.) for r η = 25, 50, thus enabling us to also obtain the l ∼ t 2 562 I. Zacharoudiou, E. M. Chapman, E. S. Boek and J. P. Crawshaw regime at early times, as the interface accelerates initially penetrating the throat. This initial acceleration of the interface for r η = 5, 25 and 50 is reflected in the increasing dynamic contact angle θ α xz , shown in figure 10(a). Rescaling length and time using (3.2) reveals that our numerical results approach the theoretical prediction of Ichikawa et al. (2004), equation (3.1), in the limit of long times as shown in figure 9(b). Here we used the equilibrium value of the contact angle θ eq = 30 • , although actually the dynamic contact angle, shown in figure 10(a), varies with time as it depends on the interfacial velocity and Ca (Cox 1986;Sheng & Zhou 1992). Using the dynamic contact angle instead would decrease slightly l * improving the agreement further. The dynamic contact angle, in the xy and xz planes, is evaluated by fitting the interface in the central region of the advancing meniscus to a circle, as shown in figure 10(c). The angle of intersection this circle makes with the side walls is θ α xy,xz . Figure 10(b) depicts the variation of cos(θ α xz ) as a function of the capillary number Ca. This is in agreement with the theoretical predicted dependency of θ α xz on Ca (Cox 1986;Sheng & Zhou 1992) cos θ α = cos θ eq − Ca ln(KL s /l s ).

(3.3)
K is a fitting constant, L s is a characteristic length scale of the system and l s is the effective slip length at the contact line. High interfacial velocities translate to high θ α xz , while as the interface slows down θ α xz approaches the equilibrium value θ eq . Fitting the results for each simulation set to (3.3) and extrapolating to Ca = 0 reveals the following θ α Ca=0 = 29.2 • , 30.3 • , 30.9 • and 29.4 • for r η = 5, 25, 50 and 500 respectively. Hence, this verifies that the dynamic contact angle tends to the correct value for the equilibrium contact angle of θ eq = 30 • . The deviation increases for the case of θ eq = 16 • and is of the order of a few degrees (∼ 5 • ). This is expected for very small or large contact angles due to the finite width of the interface and has been observed Pore-filling events in single junction micro-models for binary and ternary systems as well (Pooley et al. 2009;Semprebon, Krüger & Kusumaatmaja 2016).
We next examined the velocity of the interface front. Decreasing the viscosity ratio by decreasing the viscosity of the WP results in less viscous dissipation in the WP and hence more energy becomes available for driving the interface in the channel. Figure 11(a) shows the corresponding Ca = η w u/γ . The experimental Ca is of the order of 10 −2 (r η ∼ 50). Ichikawa et al. (2004) estimated the dimensionless velocity which, in the limit of long times, approaches u * = 1 2 cos θ eq t * . (3.5) Although they neglected variations in the advancing dynamic contact angle θ α , equation (3.5) can still be considered valid in the limit of long times as the advancing contact angle θ α approaches the equilibrium value θ eq . It is evident that the numerical results demonstrate excellent agreement with the theoretical prediction of Ichikawa et al. (2004) (figure 11b).

Junction region
Having validated the interface motion in the throat prior to the junction, we next examine the interface motion in the wider pore body. As the interface approaches the end of throat T 3 and enters the junction region, it decelerates at first as the driving capillary forces per unit area decrease and the interface adapts a concave shape in the xy-plane. This reduction in the interfacial velocity is evident in figure 11. On the other hand inertial forces can keep the interface moving in the junction, while the motion is opposed by viscous forces that can damp this forward movement. Hence, an important dimensionless number that can affect the dynamics of the interface entering the junction is the Ohnesorge number which gives the relative importance of viscous forces over inertial and capillary forces. Inertial forces (per unit volume), which scale as ∼ρu 2 /L s , are in the range 10 −6 (r η = 5)-10 −9 (r η = 500) in l.u. Viscous forces (per unit volume), which scale as ∼η w u/L 2 s , are in the range 10 −8 in l.u. for all cases as r η increases from 5 to 500. Capillary forces F cap = 2γ cos θ α (d + w t ) are of the order of 1 in l.u. for all cases.
On a second stage, the contact line in the xy-plane makes contact with the side walls, see figure 12(a). As this favours energetically a transition from a concave to a convex configuration (dashed yellow line to the solid black line configuration), the surface energy released and the interface configuration transition can lead to interfacial oscillations, which, as demonstrated in figure 12(b), are more profound with decreasing Oh. Time is normalised by the time t cont when the interface imbibes throat T 2 . Similar interfacial oscillations have been observed by Ferrari & Lunati (2013), who examined a forced imbibition situation. Here we demonstrate that these oscillations can be generated as the interface travels through a narrow throat to a wider pore body in a spontaneous imbibition scenario as well, due to the small Oh, with the mechanism behind this being the same. As the interface progresses further in the junction, the kinetic energy that was available is gradually dissipated and the (a) (3) (1) (2) (0) (4) FIGURE 13. (Colour online) Interface configuration in the xy-plane (z = d/2) in the junction for (a) r η = 5, Oh = 6 × 10 −3 and (b) r η = 50, Oh = 6 × 10 −2 .
forward movement is due to the action of capillary filling. This becomes more clear when looking at the interface configuration in the junction at a level z = d/2 for viscosity ratios r η = 5 (Oh = 6 × 10 −3 ) and r η = 50 (Oh = 6 × 10 −2 ), figure 13. Figure 14(a) shows the length travelled by the meniscus in the junction, measured along the dashed line in figure 13, for varying r η and Oh. The similar shape of the curves, with two peaks -labelled as (1) and (3) -and two troughs -points (2) and (4) -is due to the transition from a concave to convex configuration, favoured by the shape of the junction. As expected, increasing the viscosity of the wetting phase (increasing r η ) increases the time it takes for the interface to imbibe into the next downstream throat, t cont . For all situations examined here, the fluid imbibes throat T 2 first, which is achieved when the interface has progressed a length l ∼ 85 in the junction region. In other words, geometry dictates the pore-filling sequence, irrespective of the dynamics prior to the selection of the next downstream throat or the filling rules proposed by Lenormand et al. (1983). The same is observed 566 I. Zacharoudiou, E. M. Chapman, E. S. Boek and J. P. Crawshaw  figure 13(a). The wetting fluid starts imbibing the next downstream throat (T 2 ) when l ∼ 85 and t = t cont . (b) Simulation results with r η = 50 and different Oh, by varying the viscosities of both phases. Increasing fluids' viscosities increases Oh and t cont . The average Ca = η wū /γ for the interface motion in the junction is 8.1 × 10 −4 (Oh = 6 × 10 −2 ), 8.4 × 10 −4 (Oh = 2 × 10 −1 ) and 7.9 × 10 −4 (Oh = 1 × 10 0 ). The labels (0)-(4) correspond to the snapshots in figure 13(b). Inset: the moment the wetting phase starts imbibing the throat T 2 (t = t cont ) the interface retracts in the middle of the junction (reduction in l).
for the simulations reported in figure 14(b), which examines simulations with the same viscosity ratio (r η = 50) and different Oh, achieved by varying the viscosity of both phases. We note here that the capillary number is approximately the same in these simulations (Ca ∼ 8 × 10 −4 ); t cont varies significantly though as can be clearly observed, as a consequence of the different rates at which energy is dissipated in the system. Therefore, an important remark here is that Ca itself is not sufficient to characterise the fluid flow. The dimensionless numbers, relevant to the specific type of fluid flow, must be matched, for example the viscosity ratio r η and the Oh, in order to characterise the fluid flow dynamics.
Examining the imbibition process in the junction in terms of the surface energies, and given that we kept all parameters fixed except for the fluid viscosities, we note that the surface energy released, due to wetting, and used to drive the fluid-fluid interface, is the same for all runs. What changes is the rate of viscous dissipation, which dictates how much energy is left available as kinetic energy for the fluid motion. η i (i = w, nw) is the local viscosity and˙ αβ = (∂ α u β + ∂ β u α )/2 is the rate of strain tensor. For the spontaneous imbibition situation we examine here, the energy balance states (Ferrari & Lunati 2013) where E k and F are the kinetic energy and surface free energy respectively. In figure 15(a) the viscous dissipation rate is plotted as a function of time for two of the simulations reported in figure 14(b), covering the interface motion both in throat Pore-filling events in single junction micro-models for the simulations reported in figure 14 (b). This covers the fluid-fluid interface motion in the throat T 3 (t t 1 ) and in the junction region (t t 1 ). The corresponding values of viscosities are η w = 8.33 × 10 −2 , η nw = 1.67 × 10 −3 (Oh = 2 × 10 −1 ) and η w = 6.67 × 10 −1 , η nw = 1.33 × 10 −2 (Oh = 1 × 10 0 ). (b) Viscous dissipation rate versus time in scaled units.
T 3 (prior the junction) and in the junction region. In figure 15(b) viscous dissipation rate in dimensionless form Φ * = Φ/γ D hū is plotted as a function of dimensionless time t * = t/t cont . As the viscosity of both fluids increases to maintain the same r η , the total amount of energy dissipated (area under the curves, i.e. Φ dt) increases, and hence the change in kinetic energy decreases. The peak which is clearly visible for each case, corresponds to time t = t cont , when the wetting fluid starts imbibing throat T 2 , resulting in an increase in the fluid velocity.
The change in the total surface energy can be expressed as dF = γ dA int + γ ws dA ws + γ ns dA ns , (3.9) where dA int , dA ws and dA ns are the increments of the areas of the fluid-fluid, solidwetting phase fluid and solid-non-wetting phase fluid interfaces respectively and γ , γ ws , γ ns the corresponding surface tensions. Given that the total solid surface area A s tot = A ns + A ws is constant and that dA ns = −dA ws , then dF = γ (dA int − cos θ eq dA ws ). (3.10) The total surface energy is given by F = γ (A int − cos θ eq A ws ) + F 0 , where F 0 = γ ns A s tot is constant. Plotting F − F 0 as a function of dimensionless time t * = t/t cont in figure 16(a), reveals that the total surface energy decreases monotonically. It also demonstrates that the released energy is approximately the same for all simulations, as expected.
Finally, an interesting feature was observed when comparing situations with the same viscosity ratio (r η = 50), see figure 14(b). The longer times t cont observed, as viscous dissipation (and Oh) increases, allow more time to the interface to progress along the corners of the side throats. This is shown in figure 16(b), where we plot the area of the fluid-fluid interface as well as in figure 17, where we show snapshots of the interface configuration in the junction. In a complex geometry, as in porous media, swelling of these films and the consequent collapse can affect the displacement pathways. Especially in porous media and natural rock formations the 568 I. Zacharoudiou, E. M. Chapman, E. S. Boek and J. P. Crawshaw  and 15 with r η = 50 and (a) Oh = 2 × 10 −1 , (b) Oh = 1 × 10 0 . A situation with higher viscous dissipation rate (b) decreases the amount of energy converted to kinetic energy and hence the mean velocity along the x-direction. This gives more time for the development of thin wetting films progressing along the corners of the micro-channel. The snapshots correspond at the same length travelled in the x-direction and approximately at the same normalised time t * = t/t cont .
effective diameter varies continuously with length, favouring filling of the narrowest downstream throat as predicted by the filling rules of Lenormand et al. (1983).
Here, however, and examining a wide range of parameters relevant to spontaneous imbibition dynamics, LB simulations revealed that the pore-filling sequence remained the same, with the pore geometry being the major influencing factor.

Conclusions
It was observed that for both the drainage and quasi-static imbibition experiments the displacement pathways predicted by the Young-Laplace law are obeyed. LB simulations also followed these displacement predictions for drainage. In addition, the theoretical critical pressures for displacement events were in good agreement with calculated experimental events. For our spontaneous imbibition experiments, on the other hand, we found that the WP enters an adjacent throat first due to the absence of WP film growth. Lenormand (1990) suggested that the imbibing WP should enter an adjacent throat first, but no direct evidence was provided. Here we confirm this hypothesis, for the first time, directly using experiment and corresponding LB simulations. Thus, pore geometry plays a vital role as it becomes the main deciding factor in the displacement pathways. Once the critical pressure of the pore has been exceeded, all downstream throats are able to be filled. This displacement choice was observed in both our spontaneous imbibition experiments, and the corresponding LB simulations for models of the same geometry. The implications of the absence of WP films within the models have not been investigated.
Furthermore, we observe that the displacement of the meniscus in a throat and the scaling of the imbibing fluid column with time can fall in the early time flow regimes of capillary filling ((i) l ∼ t 2 , (ii) l ∼ t) prior the Lucas-Washburn regime ((iii) l ∼ sqrt (t)). An in-depth investigation of imbibition dynamics using lattice Boltzmann simulations was carried out in Zacharoudiou & Boek (2016). We emphasise here that matching the relevant dimensionless numbers is essential in correctly resolving the multiphase flow dynamics, as Ca itself is not sufficient to uniquely describe the flow. For example, we need to match the viscosity ratio and the Ohnesorge number, which for fluid flow at the pore scale is typically Oh 1. Nevertheless, for the range of parameters examined here, LB simulations revealed that the pore-filling sequence remained the same. Therefore, the main message of the current paper is that the filling order of channels connected to a junction depends primarily on the geometry of the pore body and is largely independent of the details of the dynamic meniscus shape influenced by inertial effects.
In addition to the above, in real porous media, microscopic WP films can develop under strong wetting conditions and flow through the cervices and the micro-roughness of the pore walls, whereas in our micro-fluidic devices WP films develop at the right-angled wedges. This, however, will not change the main message of the current paper. Several papers in the literature have discussed the development of thin WP films in porous media, including Vizika, Avraam & Payatakes (1994), Constantinides & Payatakes (2000), Bico & Quéré (2003). Depending on the time scales associated with the general advancement of the fluid-fluid interface (t adv ) and the time scales for thin films to develop and swell (t film ), the displacement pathways can be very different. If the fluid flow is fast, for example in the manuscript the spontaneous imbibition situation (t film > t adv ), then the displacement pathway is mainly affected by the geometry. On the other hand, if the fluid flow is very slow (providing time for thin films to develop and swell), for example in the manuscript the quasi-static imbibition situation (t film < t adv ), then the displacement pathway and the pore-filling sequence is likely to follow the filling rules proposed by Lenormand et al. (1983), based on the Young-Laplace equation.
The advancement of the invading WP fluid in a real porous medium, under strong wetting conditions, involves two distinct macroscopic fronts (Constantinides & Payatakes 2000;Bico & Quéré 2003). The primary displacement front, which saturates the medium, is due to the piston-type motion of menisci in the main pores. The secondary front is due to the precursor WP films, which propagate ahead of the primary front using the fine structures of the porous material. The distance between the two fronts depends on the capillary number Ca and the viscosity ratio and may be many times larger than the mean pore length (Constantinides & Payatakes 2000). Under certain conditions the WP films can swell and cause disconnection and entrapment of the NWP. In these situations the displacement pathways can be significantly different from situations with no WP films. The impact of this phenomenon can be significant. For example Constantinides & Payatakes (2000) report that these WP films cause a substantial increase of the residual NWP saturation. Rücker et al. (2015) report on regimes of corner flow and film swelling in real porous media, decreasing the connectivity of the NWP. Micro-fluidic devices and simulations can be used to investigate the dynamics and mechanisms involved in two phase flow at the pore scale and elucidate the role of wettability, viscosity ratio, Ca or even roughness using patterned micro-fluidic devices on the advancement of the displacement fronts.