Emergence of unstable invasion during imbibition in regular porous media

Abstract

The unstable fluid–fluid displacement patterns in porous media with rough invasion fronts and trapping of the defending phase are often observed in drainage, i.e. when the solid is non-wetting to the invading phase. On the other hand, during imbibition, compact and faceted growth is expected in regular porous media with geometrically homogeneous pore structure. Here, we report the irregular growth of invading fluid during the imbibition process in two-dimensional regular porous media. The ramified invasion morphology associated with thin fingers is reminiscent of capillary fingering. Through examining the capillary pressure signals and the type of pore-scale invasion mechanisms, the fundamental differences between faceted growth and irregular invasion are revealed. By analysing the pore-scale invasion mechanisms, a phase diagram describing the dominance of different invasion events is proposed. Through conducting systematic quasi-static radial injection simulations across a wide range of porosity and wettability, excellent agreement is observed on the transition boundary from faceted and compact displacement patterns to irregular and dendritic invasion morphologies. This is reflected by the overlap of the transition boundaries from analytical prediction, type of pore-scale invasion events, and macroscopic morphology quantified by the fractal dimension. This work provides new insights on the role of geometrical features of solid structures during multiphase flow with emphasis on the porosity and wettability. The findings could assist in guiding the design of microfluidic devices to control deterministically the multiphase flow, transport and reaction processes.

1. Introduction

The fluid–fluid displacement process is involved in various engineering applications, such as carbon geosequestration, underground hydrogen storage and fuel cells (Szulczewski et al. 2012; Andersson et al. 2016; Heinemann et al. 2021). Extensive efforts have been made to understand, predict and control the fluid transport process using simulations, experiments and theoretical models (Lenormand, Touboul & Zarcone 1988; Zhao, MacMinn & Juanes 2016; Singh et al. 2017; Rabbani et al. 2018; Guo & Aryana 2019; Gu, Liu & Wu 2021). A phase diagram of fluid–fluid displacement patterns was revealed in the seminal work by Lenormand (1990), where the three displacement regimes – including stable displacement, capillary fingering and viscous fingering – were found to be controlled by the capillary number, reflecting the relative importance of viscous force to capillary force, and the viscosity ratio of the two fluids.

Another significant progress in understanding multiphase flow in porous media is marked by the identification of impact from wettability, that is, the contact angle at the three-phase contact line measured within the invading phase (Cieplak & Robbins 1990; Trojer, Szulczewski & Juanes 2015; Jung et al. 2016; Zhao et al. 2016; Singh et al. 2017; Primkulov et al. 2018, 2019, 2021). As the solid becomes more wetting to the invading fluid, i.e. in an imbibition process, the invasion morphology shifts from dendritic patterns with many thin fingers towards a more compact region with smoother invasion front and less trapping of the defending phase. This transition was explained by the increased occurrence of the pore-scale cooperative filling events (Cieplak & Robbins 1990; Holtzman & Segre 2015), during which two menisci merge into one, consequently leading to more compact displacement patterns.

While most focus has been placed on the effects of fluid phases, the influence of topological and geometrical features of porous structures on multiphase flow received less attention despite several studies demonstrating a clear dependence of displacement patterns on the pore space geometry, especially in the quasi-static regime where capillary force controls the interface motion. It has been shown that capillary or viscous fingering can be suppressed by careful design of porous media with a particle size gradient (Rabbani et al. 2018; Lu et al. 2019). Thus understanding how the pore structure affects multiphase flow is of great importance for controlling the process of fluid transport. In regular porous media – e.g. two-dimensional micromodels composed of circular posts with narrow particle size distribution arranged on a triangular or square lattice – both simulations and experiments have shown compact and faceted patterns (Cieplak & Robbins 1990; Martys, Robbins & Cieplak 1991; Holtzman 2016; Wang, Pereira & Gan 2021). This phenomenon was attributed to the dominance of the underlying pore geometry (regular lattice structure), which results in the symmetrical growth of invasion patterns. However, when heterogeneity is introduced to particle shape (Wang et al. 2021), size (Lenormand 1990; Holtzman 2016; Hu et al. 2019) or position (Wang et al. 2019), irregular displacement patterns with rough invasion fronts are observed. These observations can be explained by the amplified unevenness in the local critical capillary pressures due to heterogeneity in pore structure, leading to ramified advancement of menisci even in the wetting condition (imbibition) (Wang et al. 2019). One limitation of the aforementioned works, as pointed out by Golmohammadi et al. (2021), is that the porosity of porous media is limited to very high values (0.6–0.9) with sparsely packed particles, whereas the multiphase behaviour in porous media with lower porosity remains unexplored.

Clearly, there is a need for thorough examination of the role of microstructure during multiphase flow in porous media. In this work, we study the imbibition process in regular porous media across a wide range of porosity in the quasi-static regime. As expected, due to homogeneity of the porous structure and the wetting condition during imbibition that favours cooperative pore filling events, the invasion morphology is compact and faceted (Lenormand 1990; Holtzman 2016; Gu et al. 2021; Wang et al. 2021). Interestingly, when the porosity of porous medium decreases below a critical value, the emergence of unstable invasion associated with rough invasion fronts and occurrence of trapping events is observed, which is reminiscent of invasion percolation (or capillary fingering in drainage condition). By analysing the pore-scale instability events, we show that the macroscopic transition of multiphase flow from a compact and faceted displacement towards an irregular and unsymmetrical growth is linked to the shift in the type of pore filling events. Based on this, a phase diagram is proposed to predict different regimes of multiphase displacement patterns in regular porous media. Results from systematic fluid–fluid displacement simulations with different contact angles and porous structures agree well with the proposed phase diagram.

2. Method

A pore-resolved interface tracking algorithm is used for simulating multiphase flow processes in the capillary-dominated regime (Wang et al. 2021, 2022). The capillary-dominated regime (in the limit of vanishing capillary number), or the quasi-static regime, corresponds to displacement processes where the viscous effect (regardless of viscosity ratio) is negligible compared with the interfacial tension. The algorithm incorporates different types of pore-scale instabilities during fluid invasion, including burst, touch and overlap events (figures 1a–c) (Cieplak & Robbins 1990). The solid surface is discretized into computational nodes such that the sharp edge pinning effect can be taken into account (Gibbs 1961; Oliver, Huh & Mason 1977). The capillary pressures associated with different meniscus advancement events are determined based on the two-dimensional Young–Laplace equation, i.e. assuming that the contribution from the out-of-plane curvature can be neglected: $P_c=\gamma /R$, with $\gamma$ the interfacial tension and $R$ the in-plane meniscus radius. At each computation step, the meniscus with the minimum capillary resistance advances. The algorithm tracks the interface motion during the fluid–fluid displacement process, and the simulation stops when the invading fluid reaches the outlet. The evolution of capillary pressure can also be obtained during the movement of menisci in pore space. Validation of the numerical method in a previous study has indicated that the invasion mechanisms at pore scale are captured accurately, and the macroscopic invasion morphology of fluid–fluid displacement in disordered porous media can be reproduced. The key advantage of the algorithm is its applicability to multiphase flow in porous media with arbitrary pore structures. However, it is currently limited to capillary-dominated processes, and further efforts are needed to incorporate viscous effects as well as extending the algorithm to three dimensions. A detailed description of the algorithm can be found in the previous study (Wang et al. 2021, 2022) and is omitted here for brevity.

Figure 1. Fluid displacement in regular porous media composed of circles packed on a triangular lattice. (a) Schematic of a burst invasion event, where the meniscus reaches the local maximum capillary resistance pressure and becomes unstable. (b) A touch event, where the meniscus comes into contact with a particle, forming two new menisci. (c) An overlap event, where two menisci intersect and merge into one meniscus. (d,e) Invasion morphologies at percolation (when the invading phase reaches the boundary) during radial injection simulations for porous media with porosity $\phi =0.65$ and 0.25, respectively. The invading phase saturations are 0.275 and 0.395 for (d,e), respectively. Both cases have intrinsic contact angle $\theta _0=60^\circ$. The insets show magnified regions of the porous media for each case.

Two types of regular lattice structures, triangular and square lattices, are considered in this work. For both cases, the porous media are composed of about 2000 circular posts of equal size arranged at corresponding lattice sites. Each circle with radius unity is discretized into 100 computational nodes, i.e. represented by a 100-edge regular polygon, after which a random rotation about the particle centre by an angle ranging between $0^\circ$ and $360^\circ$ is applied. This manually introduced rotation, which can be regarded as a small perturbation in geometry, prevents the exact equality of capillary pressure of different menisci in a perfectly crystalline lattice with a single post radius. Mesh independence studies were conducted to ensure that the mesh size is sufficiently small. Radial injection simulations are conducted for porous media with different porosities, and the simulation stops at percolation, i.e. when the invading phase reaches the perimeter (boundary particles) of the domain. The explored intrinsic contact angles for imbibition processes range from $45^\circ$ to $90^\circ$. Since the current algorithm does not capture the phenomenon of corner flow in the extreme wetting condition (Zhao et al. 2016), cases with contact angle $\theta _0<45^\circ$ are not considered in the present work.

3. Results and discussion

3.1. Symmetry breaking in pattern formation

Figures 1(d,e) show the final invasion morphologies of invading phase (blue) with $\theta _0=60^\circ$ and porosity $\phi =\{0.65,0.25\}$, respectively, with insets showing the corresponding particle packing structure of a magnified sub-region. At $\phi =0.65$, the pattern exhibits a hexagonal shape due to the dominance of the underlying triangular lattice on which the particles are arranged, consistent with past experiments and numerical simulations (Lenormand 1990; Holtzman 2016; Wang et al. 2021). For $\phi =0.25$, however, a ramified shape with several clusters of trapped defending phase is seen, manifesting the transition towards the regime of capillary fingering described by the theory of invasion percolation (Wilkinson & Willemsen 1983; Lenormand & Zarcone 1985). Extra simulations for several different realizations of porous media with the same key parameters, i.e. $\theta _0$ and $\phi$, as well as different mesh sizes (polygon edge number 50, 100 and 200), are conducted, and the same qualitative results are observed.

To illustrate the differences observed in figures 1(d,e) from another perspective, we plot the capillary pressure as well as corresponding pore-scale instability events during the displacement process in figure 2. We define the dimensionless capillary pressure as $P_c^*=a/R$, with $R$ the meniscus radius and $a$ the particle centre-to-centre spacing. In figure 2(a), the $P_c^*$ evolution during a sub-period of the entire displacement process for the porosity $\phi =0.65$ case shows a periodic behaviour where the touch events associated with large $P_c^*$ are followed by multiple alternating overlap and touch events, before the next cycle starts. Note that the data shown in figure 2 are representative, i.e. the periodicity persists outside the presented time window, and the length of each cycle (number of these alternating events) increases as the invasion area grows larger. By visualizing the invasion progress (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2022.336), it is found that this periodic behaviour is associated with consecutive filling of pores along an edge of the hexagon; i.e. starting from a hexagon, once a new pore is filled, the invasion stops along all other edges, until the current edge is filled completely. The distribution of $P_c^*$ values reveals that there are in total three peaks. The peak with large $P_c^*$ is associated with relatively small number of touch events. Note that the fluctuations of $P_c^*$ within each peak are attributed to the geometrical noise introduced by the random rotation of polygons. The non-overlap of peaks indicates that the geometrical perturbation is small enough, and the mesh size (circular posts represented by 100-edge polygons) is sufficiently fine, such that the local invasion mechanisms in regular lattice with monodisperse particles are preserved. As expected, it is found that the fluctuation decreases as the number of mesh points increases.

Figure 2. Capillary pressures and invasion events during the displacement process. (a) Evolution of the dimensionless capillary pressure $P_c^*$ with overlap and touch events marked by red circles and blue squares, respectively. The result corresponds to the case in figure 1(d), i.e. $\theta _0=60^\circ$, $\phi =0.65^\circ$. Only a sub-period of the total capillary signal is shown for clarity. The plot on the right shows the distribution of $P_c^*$ values corresponding to overlap (red) and touch (blue) events, respectively. The number of events during the entire simulation is shown in the legend. (b) Data for the case with $\theta _0=60$, $\phi =0.25$ (case in figure 1e).

For $\phi =0.25$, figure 2(b) shows that the overlap events occur intermittently between touch events that are all associated with high capillary pressure. Also, in contrast to the previous case, the distribution of $P_c^*$ contains only two peaks, a clear indication of change in the capillary pressure signature. The invasion morphology is found to evolve according to a rather chaotic manner with a dendritic pattern, reminiscent of the process of invasion percolation (see supplementary movie 2). The asymmetrical advancement of invading fluid, i.e. symmetry breaking in the pattern formation, is intuitively unexpected since the porous medium exhibits a high degree of homogeneity, and the contact angle $\theta _0=60^\circ$ falls in a range (imbibition process) where the cooperative pore filling events are favoured.

3.2. Analysis of critical pore-scale events

To explain the observed regime transition, we analyse the pore-scale invasion mechanisms for regular porous media at different contact angles. As a meniscus moves forward between a pair of circular pillars, its radius can be calculated as (Purcell 1950)

(3.1)\begin{equation} R=\frac{a-2r\cos{\alpha}}{2\sin{(\alpha+\theta_0-90^\circ})}, \end{equation}

with $a$ the particle centre-to-centre spacing, $r$ the particles radius, $\alpha$ the angle associated with the meniscus position (inset in figure 3a). Figure 3(a) shows the value of $P_c^*$ at different $\alpha$ and $\theta _0$ with ${r}/{a}=\frac {5}{16}$, corresponding to $\phi =0.65$ (case in figures 1(d) and 2(a)). As the meniscus advances (from left to right, i.e. increasing $\alpha$), it may experience pore-scale instabilities including burst, touch or overlap (Cieplak & Robbins 1990). The $P_c^*$ value associated with burst is the maximum value marked by triangles in figure 3(a), which also gives the corresponding value of $\alpha _{burst}$ (the value of $\alpha$ at which burst occurs). The configuration of the meniscus shape corresponds to an arc whose centre is located at the midpoint of the segment that connects the particle centres. The touch event occurs when the meniscus comes into contact with the particle in the front, which satisfies geometrically (see Appendix A for the derivation)

(3.2)\begin{equation} \left[1-\cos{(\alpha_{touch}+\theta_0-90^\circ)}\right]R+r\sin{\alpha_{touch}}=\tfrac{\sqrt{3}}{2}a-r. \end{equation}

After replacing the meniscus radius $R$ according to (3.1), equation (3.2) involves only one variable, $\alpha _{touch}$, and can be solved numerically. The determined $\alpha _{touch}$ and its corresponding dimensionless capillary pressure $P_c^*=a/R$ are marked by squares in figure 3(a). The value of $\alpha$ for an overlap event can be obtained from (Appendix A)

(3.3)\begin{equation} \alpha_{overlap}=\max[\theta_{c},2\theta_{c}-\min(\alpha_{touch},\alpha_{burst})], \end{equation}

with $\theta _c=60^\circ$ the characteristic angle of porous media with triangular lattice (black dashed line in figure 3a) (Holtzman & Segre 2015). Equation (3.3) implies that $\alpha _{overlap}$ depends not only on the intrinsic lattice structure ($\theta _{c}$), but also on the position of its neighbour meniscus before reaching an instability event.

Figure 3. Determination of critical pore-scale events. (a) The dimensionless capillary pressure $P_c^*$ as a meniscus advances in the porous medium. The location of the meniscus is indicated by $\alpha$ in the inset. The porous medium has porosity $\phi =0.65$ (corresponding to $r/a=5/16$). The burst, touch and overlap events are marked by triangles, squares and circles, respectively. A filled marker indicates the first instability event (critical invasion event) that a meniscus experiences. Plots of $P_c^*$ before and after these critical events are drawn in solid and dotted lines, respectively. The black dashed line marks the characteristic angle for overlap, which is $60^\circ$ for circles placed on a triangular lattice. (b) Plots of $P_c^*$ and invasion events for the case with smaller porosity $\phi =0.25$. (c) Phase diagram showing the type of critical invasion events for different porosity and wettability. The red dotted line and red dashed line correspond to results shown in (a,b), respectively.

As a meniscus moves forward during the invasion process, it is the first critical instability, not necessarily the one associated with the smallest $P_c^*$, that determines whether a local pore will be invaded. The $P_c^*$ associated with the first instabilities are represented in figure 3(a) by filled markers, and the capillary pressures determined from (3.1) before and after touch or overlap events are represented by solid and dashed lines, respectively. Thus the dashed line represents the theoretical but not realistic capillary pressure that a meniscus can never experience due to the occurrence of merging/splitting events. Figure 3(a) shows that as the wetting condition changes from drainage ($\theta _0>90^\circ$) to imbibition ($\theta _0<90^\circ$), the critical pore-scale invasion mechanisms shift from burst to overlap. For smaller porosity $\phi =0.25$ (figure 3b), it can be seen that $P_c^*$ is amplified in both imbibition and drainage due to narrower throats. Also, smaller values of $\alpha$ associated with the maximum $P_c^*$ (burst events) are observed, indicating that the overlap events become less favourable due to the decrease of porosity. A consequent effect is the emergence of touch events (solid squares) in the weak-imbibition regime, which are absent in sparser media (figure 3a). By sweeping the $\theta _0$–$\phi$ space, a phase diagram characterizing the critical type of pore-scale instability can be established (figure 3c). The red dotted line and red dashed lines correspond to figures 3(a) and 3(b), respectively. It has been identified that the overlap event, or the cooperative pore filling event (Cieplak & Robbins 1988; Holtzman & Segre 2015), smooths the invasion front and makes the invasion morphology compact, whereas the occurrence of either burst or touch events destabilizes the displacement. Clearly, crossing of regimes is possible by varying either $\theta _0$ or $\phi$. Specifically, a smaller $\theta _0$ and greater $\phi$ favour the overlap events, whereas greater $\theta _0$ and smaller $\phi$ promote touch or burst events.

3.3. Phase diagram

The main results of this work are shown in figure 4, which compares the systematic numerical simulations with the analytical prediction in the imbibition regime. The final invasion morphologies (invading fluid in blue colour) of radial injection simulations at percolation are shown in figure 4(a) for the intrinsic contact angles $\theta _0\in [45^\circ,90^\circ ]$ and $\phi \in [0.2,0.85]$. Here, particles that are fully submerged in the invading fluid are also coloured in blue. Qualitatively, as $\theta _0$ increases and $\phi$ decreases, a clear transition from compact hexagonal patterns towards dendritic shapes with rough invasion fronts is observed. The blue and brown boxes correspond to the cases in figures 1(d) and 1(e), respectively. To quantify the displacement patterns, in figure 4(b), we present (i) the proportion of overlap events $P_{overlap}$ during the entire displacement process, (ii) the fractal dimension $D_{f}$ of final invasion morphologies using the box-counting method (Mandelbrot 1983), and (iii) the boundary of overlap obtained analytically in figure 3(c) and the corresponding adjusted line for consideration of the sharp edge pinning effect (Gibbs 1961; Oliver et al. 1977; Wang et al. 2021) by shifting the original curve towards the left by $3.6^\circ$. This is because each particle is represented by a 100-edge polygon that has a local corner angle of $176.4^\circ$, leading to an effective contact angle $\theta =\theta _0+(180^\circ -176.4^\circ )=\theta _0+3.6^\circ$ during simulations (Appendix A).

The three metrics described in (i)–(iii), respectively, represent the pore-scale, macro-scale and analytical characterization of the displacement processes. Below the adjusted analytical curve, the overlap events account for more than 40 % of the total pore invasion events, and the corresponding fractal dimension remains close to 2 (note the adjusted colourbar for $D_f$, where the brown colour represents values of around 1.98). Above the analytical curve, both $P_{overlap}$ and $D_{f}$ experience a sharp decrease. Overall, figure 4(b) demonstrates good consistency among these multi-scale descriptions. The adjusted analytical curve (black-dashed line) is mapped onto figure 4(a), which again demonstrates the accurate prediction of the phase transition.

Figure 4. Phase diagram for displacement patterns during imbibition in regular porous media. (a) Invasion morphologies of radial injection simulations at percolation for different wettability and porosity. The black dashed line is mapped from the black dashed line in (b). The blue and brown boxes correspond to the cases in figure 1(d) (also figure 2a) and figure 1(e) (also figure 2b), respectively. (b) Proportion of overlap events $P_{overlap}$ contour. The colours of dots represent the measured fractal dimension $D_f$ of the invasion morphology using the box-counting method. The black dotted line marks the boundary of overlap events (figure 3c). The black dashed line is obtained by shifting the black dotted line by $3.6^\circ$ to the left for consideration of the sharp edge pinning effect due to the finite number of mesh points.

3.4. Generalization to square lattice

We then generalize our results to porous media with particles packed on a square lattice, which is another common regular configuration. Analytically, the calculation for a burst event remains the same, and (3.1) remains valid. For a touch event, the value of $\alpha _{touch}$ follows (Appendix A) from

(3.4)\begin{equation} \left[R\cos{(\alpha_{touch}+\theta_0-90^\circ)}+a-r\sin{\alpha_{touch}}\right]^2+ \left(\frac{a}{2}\right)^2=(R+r)^2, \end{equation}

again with $R$ expressed in (3.1). Equation (3.3) remains applicable for overlap events, but the characteristic angle for the square lattice becomes $\theta _{c}=45^\circ$ (black dashed line in figure 5a). The $P_c^*$ evolution as a function of meniscus position $\alpha$ is shown in figure 5(a) for different intrinsic contact angles. The $\theta _0$–$\phi$ phase diagram is shown in figure 5(b). Compared with figure 3(c) for a triangular lattice, one can notice the absence of touch events. This is expected intuitively since the parallel alignment of particles (inset in figure 5a) suppresses the occurrence of a touch event without significant advancement of meniscus position ($\alpha$), whereas a staggered particle arrangement would favour comparatively the touch events. Also, a slight expansion of the overlap region (from $99^\circ$ in figure 3(c) to $110^\circ$ in figure 5(b) at $\phi =0.85$) is observed. This is attributed to a smaller $\theta _c$, resulting in more favoured overlap events.

To check the validity of the analysis and also preserve brevity, only six simulations with $\phi =\{0.3, 0.5, 0.8\}$ and $\theta _0=\{60^\circ,90^\circ \}$ near the phase transition are conducted (grey dots in figure 5b), and the corresponding displacement patterns at percolation are shown in the adjacent insets. In the overlap regime, the patterns have rectangular shapes, revealing the dominance of the underlying lattice structure. In the burst regime, all invasion morphologies exhibit irregular patterns with thin fingers and trapping of the defending phase. These observations validate the derived phase diagram, showing the generality of the analysis.

Figure 5. Generalization to a square lattice. (a) Plots of $P_c^*$ and invasion events for a porous medium with a square lattice configuration. Symbol meanings are the same as in figure 3(a). The characteristic angle of overlap for the square lattice is $45^\circ$, as denoted by the black dashed line. (b) Invasion events phase diagram for the square lattice. The red dashed line corresponds to the path of (a). Insets show the invasion morphologies at percolation for $\phi$ and $\theta _0$, represented by corresponding dots.

3.5. Discussions

Several comments can be made based on the insights provided in this work. Recently, Gu et al. (2021) proposed a mathematical model to investigate the imbibition processes in a dual-permeability pore network, where the circular particles are arranged on a square lattice with porosity $\phi =0.55$ and $\theta _0=30^\circ$. This configuration, according to figure 5(b), falls into the overlap regime. As expected, Gu et al. (2021) observed compact growth of the invading fluid. They further claimed the independence of porosity and contact angle of their results by conducting extra simulations at $\phi =0.62$ and $\theta _0=\{15^\circ,45^\circ \}$. However, these settings again imply an overlap-dominated regime according to the phase diagram proposed in the present work. We hypothesize that the growth of invading fluid during the imbibition process could be different qualitatively in the dual-permeability network if their simulations were conducted at a much lower porosity.

In the current work, regular media with mono-sized particles are considered. If geometrical heterogeneity is introduced – for example, by considering either a distribution of particle sizes or an insufficient number of mesh points for describing the particles (or limited resolution in microfluidic devices manufacture process for experiments) – then a wider transition across the phase boundary as in figure 4(b) can be expected. Primkulov et al. (2018) studied the quasi-static fluid displacement processes in a disordered porous media, and demonstrated that an increase in posts spacing, i.e. increase in porosity, widens the range of contact angles for cooperative filling events, making the cooperative pore filling more favourable, leading to more compact displacement. One possible way to evaluate quantitatively the ‘irregularity tolerance’ of porous media for achieving compact growth could be derived by analysing the stability of the occurrence of critical invasion events, i.e. the distance (difference in $\alpha$ in, e.g. figure 3a) between the first and second instability events. From figures 3(a) and 3(b), it can be shown that the ‘irregularity tolerance’ is greater for smaller $\theta _0$ and greater $\phi$, which corresponds to moving towards the bottom left corner in the $\theta _0$–$\phi$ phase diagram. This may explain the observed hexagonal invasion morphology during imbibition even with 52 % variation in particle size (Holtzman 2016).

4. Conclusions

We have investigated quasi-static multiphase flow during imbibition processes in two-dimensional regular porous media. It is found that for relatively sparse porous media (high porosity), the displacement process is similar to crystal growth, where the invading fluid remains compact and faceted with its shape governed by the underlying lattice structure. However, when the porosity decreases below a critical value, the invasion morphologies shift towards ramified patterns associated with thin fingers and trapping of the defending fluids, which is reminiscent of capillary fingering during drainage processes. Through analysing different invasion mechanisms at pore scale, the transition from faceted to irregular growth, as the porosity decreases, is explained by less favoured overlap events. The analysis across a wide range of porosity and wetting conditions leads to a phase diagram describing the dominance of different types of invasion events, which can be used to predict the macroscopic morphological transition. By conducting systematic imbibition simulations on regular porous media of different lattice structures, i.e. triangular and square lattices, excellent consistency is observed among pore-scale events, macroscopic invasion morphologies and analytical predictions.

This study provides new insights on the role of geometrical features of solid structures during the fluid–fluid displacement process, with emphasis on the porosity and wettability of regular porous media. The findings would assist in guiding the design of microfluidic devices to control deterministically the multiphase flow, transport and reaction processes.