1. Introduction
Multiphase flow in porous media has attracted significant research interest due to its critical role in various engineering applications, including soil and underground water remediation (Karadimitriou et al. Reference Karadimitriou, Joekar-Niasar, Babaei and Shore2016), subsurface carbon dioxide and hydrogen storage (Krevor et al. Reference Krevor, de Coninck, Gasda, Ghaleigh, de Gooyert, Hajibeygi, Juanes, Neufeld, Roberts and Swennenhuis2023) and hydrocarbon recovery (Blunt Reference Blunt2017). Although extensive studies have investigated the role of wettability and its significant impacts on multiphase flow pattern formation (Singh et al. Reference Singh, Jung, Brinkmann and Seemann2019; Dalton, Crandall & Goodman Reference Dalton, Crandall and Goodman2020), less attention has been paid to the complex interaction between wettability and geometric confinement in porous media with different dimensions. Pioneering pore-scale insights into the wettability effects on displacement mechanisms were derived primarily from (quasi-) two-dimensional (2-D) microfluidic chips under piston-type injection, which largely presented a monotonic wettability effect. The mechanism lies in that, as the system transitions from drainage to imbibition (where the invading fluid contact angle
$\theta \lt 90^\circ$
), compact growth of the fluid–fluid interface is expected due to cooperative pore filling (Stokes et al. Reference Stokes, Weitz, Gollub, Dougherty, Robbins, Chaikin and Lindsay1986; Cieplak & Robbins Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990; Holtzman & Segre Reference Holtzman and Segre2015; Jung et al. Reference Jung, Brinkmann, Seemann, Hiller, Sanchez de La Lama and Herminghaus2016; Singh et al. Reference Singh, Scholl, Brinkmann, Michiel, Scheel, Herminghaus and Seemann2017). These pore-scale events generally support the monotonic wettability rule that reducing the contact angle of the invading fluid can change the incomplete displacement to a complete one, thereby improving the displacement efficiency. With the development of experimental technology, more complex porous media have been designed and manufactured, such as chromatography columns, glass beads, membranes and structured microfluidic chips, and the effect of wettability on displacement efficiency becomes a debate: some studies report a monotonic relationship while others observe a non-monotonic one (Jung et al. Reference Jung, Brinkmann, Seemann, Hiller, Sanchez de La Lama and Herminghaus2016; Zhao, MacMinn & Juanes Reference Zhao, MacMinn and Juanes2016; Singh et al. Reference Singh, Scholl, Brinkmann, Michiel, Scheel, Herminghaus and Seemann2017; Bakhshian et al. Reference Bakhshian, Rabbani, Hosseini and Shokri2020; Zulfiqar et al. Reference Zulfiqar, Vogel, Ding, Golmohammadi, Küchler, Reuter and Geistlinger2020; Lei et al. Reference Lei, Lu, Gong and Wang2023a
, Reference Lei, Lu, Liu and Wang2022; Wang et al. Reference Wang, Pereira, Sauret and Gan2023). Currently, understanding and manipulation of the multiphase displacement dependent on wettability in porous media remains an open and underexplored topic.
Previous studies have attributed the non-monotonic wettability effect mainly to the heterogeneity of geometries or flow conditions in multiphase displacements in porous media. For example, a non-monotonic wettability effect on multiphase displacement was reported in homogeneous 2-D microfluidic porous medium with asymmetric inlet–outlet conditions under viscously unfavourable fluid–fluid displacement (viscosity ratio
$M={\mu }_i/{\mu }_d \ll 1$
, where
$\mu _i$
and
$\mu _d$
are the dynamic viscosity of the invading and defending fluids, respectively) and injection from the centre of the microfluidic disk (Trojer, Szulczewski & Juanes Reference Trojer, Szulczewski and Juanes2015; Zhao et al. Reference Zhao, MacMinn and Juanes2016). The viscous instability and the point-injection condition produced the flow heterogeneity in the homogeneous porous material. Even under the stable displacement conditions (
$M \gt 1$
), the non-monotonic wettability behaviour was observed by shifting from the traditional piston-type injection to a preferential flow boundary condition (Lei et al. Reference Lei, Lu, Gong and Wang2023a
, Reference Lei, Lu, Liu and Wang2022). However, by using porous media with three-dimensional (3-D) geometries, such as glass bead models, the non-monotonic effect of wettability on multiphase displacement was also obtained, even without any heterogeneity of the geometries or flow conditions, namely under viscous-stable conditions and piston-type injection in homogeneous porous media (Zulfiqar et al. Reference Zulfiqar, Vogel, Ding, Golmohammadi, Küchler, Reuter and Geistlinger2020; Geistlinger et al. Reference Geistlinger, Golmohammadi, Zulfiqar, Schlueter, Segre and Holtzman2024). These pore-scale studies indicate that non-monotonic wettability effects are influenced not only by flow or geometric heterogeneity but also by geometric confinement when the dimension of the model system changes. Although efforts have been made to mimic the 3-D pore geometry in 2-D microfluidic chips, which is continuing to advance (Yun et al. Reference Yun, Ross, Roman and Kovscek2017; Lei et al. Reference Lei, Lu, Gong and Wang2023b
, Reference Lei, Gong, Lu and Wang2024), the crucial impact of the 3-D pore geometry on the non-monotonic wettability effect in porous media remains unclear.
To better understand the pore-scale events observed in the 2-D models which may contribute significantly to the real 3-D flows and displacements in porous media, we introduce a hybrid model, 2.5-dimensional (2.5-D) porous structures on a chip, in which the pore space contains the 3-D geometric features (e.g. independent aspect ratios and weak depth confinement), but the topology of network remains planar as in the 2-D model for visualisation, as shown in figure 1. By combining microfluidic experiments and a network model, we demonstrate that a monotonic wettability effect is observed in 2-D porous media, while the 3-D pore-geometry confinements in 2.5-D and 3-D porous structures can trigger a non-monotonic wettability effect. Direct numerical simulations revealed that the
$Ca$
–
$\theta$
phase diagram becomes more complex due to 3-D geometry. The capillary number is defined as
$Ca={\mu _i}U/\gamma$
, where
$U$
is the characteristic flow velocity and
$\gamma$
is the interfacial tension. The transitions of these pore-scale events will lead to different displacement patterns and different wettability effects. Our findings show how the 3-D pore geometry influences optimal wetting conditions to maximise displacement efficiency, which can be favourable to control displacement in natural porous media (e.g. rock), as well as offering geometric design principles for artificial porous media (e.g. microchips).
Comparison of 2.5-D porous media with 2-D and 3-D porous media in terms of (a) network topology, (b) pore-geometry features and (c) pore-size distributions.
2. Methods
Figure 1 indicates the 2-D, 2.5-D and 3-D porous media and summarises their differences, which will be used in our pore-scale studies. The 3-D porous medium is generated from X-ray micro-computed tomography (CT) images of a rock sample from the Changqing sandstone formation in Shaanxi, China, with a resolution of 0.8
$\unicode{x03BC}$
m pixel−1. To accurately reproduce the statistical information inherent in 3-D rocks, such as the pore-size distribution and the hierarchical pore structure, the microfluidic porous geometry will be generated using the random generation algorithm (Lei et al. Reference Lei, Liu, Xie, Yang, Wu and Wang2020). Using our recent innovative fabrication technique (Lei et al. Reference Lei, Lu, Gong and Wang2023b
, Reference Lei, Gong, Lu and Wang2024), the features of the 3-D pore geometry can be incorporated into the 2.5-D microfluidic porous media. Due to the small pore structures and complex porous geometries, these microfluidic chips will be fabricated based on silicon material. Detailed fabrication methods are introduced in supporting information, Text S1 and figure S1.
A microfluidic visualisation platform was meticulously designed to conduct multiphase displacement experiments in both 2-D and 2.5-D porous media. Detailed information on microfluidic visualisation systems is introduced in the supporting information, Text S1 and figure S1. For the strong imbibition condition (
$\theta = 27^\circ \pm 3^\circ$
), water was used as the invading fluid, air as the defending fluid, and the solid surface was silica, achieved through thermal oxidation on a silicon wafer (Chomsurin & Werth Reference Chomsurin and Werth2003). For the weak imbibition condition (
$\theta = 62^\circ \pm 3^\circ$
), water again acted as the invading fluid, with decane as the defending fluid, and the solid surface remained silica through thermal oxidation on the silicon wafer. In the drainage condition
$\theta = 119^\circ \pm 5^\circ$
, water served as the invading fluid, decane as the defending fluid, and the solid surface was the original silicon wafer, which had been stored in nitrogen for over a month. The dynamic viscosities of water, air and decane are 1.0 cp, 0.015 cp and 0.84 cp, respectively. Viscous instability can be neglected in this study due to the viscosity ratio
$M \gt 1$
(Lenormand, Touboul & Zarcone Reference Lenormand, Touboul and Zarcone1988). The interfacial tension
$\gamma$
for the water–air interface is 72 mN/m, and 49 mN/m for the water–decane interface. Image analysis is introduced in supporting information, Text S2 and figure S2.
Multiphase flow processes in porous media can be captured by solving Navier–Stokes (N-S) equations in the pore space and employing the volume-of-fluid (VOF) method to track the evolution of the fluid–fluid interface (Ferrari & Lunati Reference Ferrari and Lunati2013). In our simulations, the direct solution of these governing equations was carried out using the open-source software, OpenFOAM. The interFoam solver is adopted for multiphase flow simulations. Gravity was neglected in this study due to the small Bond number (
$Bd \ll 1$
). Details about the direct numerical simulations and the corresponding benchmarks can be found in the supporting information, Text S3 and figure S3. Given the significant computational costs associated with these simulations, we focused on a local region of the microfluidic chip (2 mm
$\times$
1.5 mm) as the computational domain, which meets the criteria of a representative elementary volume based on porosity evaluations. Even with this localised porous structure, the computational demands remain substantial with approximately 5.5 million grid blocks.
An intuitive and computationally inexpensive alternative for simulating porous flows is the network models, where the porous geometry is simplified into interconnected pores and throats (Joekar-Niasar & Hassanizadeh Reference Joekar-Niasar and Hassanizadeh2012). Compared with a quasi-steady network model determined by capillary forces only, the dynamic network models (DNMs) include interfacial events controlled by both viscous and capillary forces. Once interfacial events are updated based on the local pressure and velocity through their constitutive models, the overall pressure and velocity fields can be determined using the conservation equations – Kirchhoff’s equations. The DNM has been successfully established by incorporating pore-scale invasion events, such as main meniscus flow, corner flow, snap-off, by-pass, and viscous-capillary coupling effect on hydraulic conductance, in our previous studies (Gong et al. Reference Gong, Liu, Lei, Ju and Wang2023; Lei et al. Reference Lei, Lu, Gong and Wang2023a , Reference Lei, Lu, Gong and Wangb , Reference Lei, Gong, Lu and Wang2024). Further details and the corresponding benchmarks of the DNM are provided in supporting information, Text S4 and figure S3.
The transition from monotonic to non-monotonic wettability effect across 2-D to 3-D porous media. (a) The DNM simulation results about representative multiphase distribution in the 2-D, 2.5-D and 3-D porous media from strong imbibition to drainage conditions. (b) Flow rate distribution in 2-D, 2.5-D and 3-D porous media. (c,d) Variation in the invading fluid saturation at breakthrough versus the contact angle
$\theta$
in (c) 2-D, (d) 2.5-D and 3-D porous media. The inset images depict the displacement patterns of microfluidic experiments in (c) 2-D and (d) 2.5-D porous media under
$\theta =30^\circ$
,
$60^\circ$
and
$120^\circ$
. The y-axis error bar is generated from the three repeated displacement experiments.
3. Results and discussion
3.1. Wettability effect transition from 2-D to 3-D porous media
In this study, we assume a uniform distribution of wettability throughout the porous materials. The corresponding capillary number is
$Ca = 2.6\times {10^{-6}}$
. A series of DNM simulations are performed in 2-D, 2.5-D and 3-D porous media with a wide range of contact angles (
$20^\circ \leqslant \theta \leqslant 160^\circ$
), as shown in figure 2. Figure 2(a) presents the representative phase distribution of the invading fluid at the breakthrough time in 2-D, 2.5-D and 3-D porous media with three different contact angles (
$\theta = 20^\circ , 60^\circ$
and
$120^\circ$
). The phase diagram demonstrates the dependence of the morphology of invasion patterns on the wettability condition and the geometric confinement of porous media. For 2-D porous media, a decrease in
$\theta$
yields a more stable, compact displacement with more uniform and intertwined network invading patterns (figure 2
a). However, a reversal in the wettability effect on displacement patterns is observed for 2.5-D and 3-D porous media. Decreasing
$\theta$
results in a more stable and compact displacement of the defending fluid until a critical wetting transition, after which the displacement becomes significantly unstable and incomplete.
The invading saturation variation
$S_i$
at the breakthrough stage was defined as the displacement efficiency
$E_d$
(figures 2
c and 2
d). The inset images and star symbols in figures 2(c) and 2(d) show the displacement patterns and
$E_d$
obtained by the microfluidic experiments in 2-D and 2.5-D porous media. The dynamic network modelling can be validated by comparing displacement patterns and
$E_d$
(figures 2
c and 2
d) with the microfluidic experiments. Compared with traditional DNMs, our improved DNM can more accurately reflect the microscopic fluid dynamics and macroscopic multiphase flow pattern, Text S5 and figures S4 and S5. A closer inspection of the results further demonstrates that the invading fluid saturation profile of 2-D to 3-D porous media is strongly linked to the displacement patterns. The results for the 2-D porous media (figure 2
c) indicate a monotonic behaviour for
$E_d$
with
$\theta$
. The
$E_d$
continuously increases as the wetting condition changes from drainage to imbibition. This observation aligns with traditional monotonic wettability rules (Stokes et al. Reference Stokes, Weitz, Gollub, Dougherty, Robbins, Chaikin and Lindsay1986; Cieplak & Robbins Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990; Holtzman & Segre Reference Holtzman and Segre2015; Jung et al. Reference Jung, Brinkmann, Seemann, Hiller, Sanchez de La Lama and Herminghaus2016; Singh et al. Reference Singh, Scholl, Brinkmann, Michiel, Scheel, Herminghaus and Seemann2017). Interestingly, in contrast to the 2-D porous media,
$E_d$
is a non-monotonic function of
$\theta$
in 2.5-D and 3-D porous media, indicating maximum displacement efficiency under weak imbibition conditions. The wettability effects in 2.5-D and 3-D porous media are similar, which highlights that the 3-D pore geometry, rather than the network topology, plays a crucial role in the non-monotonic wettability effect (figure 2
d). According to the flow field distribution, 2.5-D and 3-D porous media contain stronger flow heterogeneity compared with the 2-D porous media (figure 2
b). This observation aligns with our recent work that preferential flow can induce the non-monotonic wettability effects (Lei et al. Reference Lei, Lu, Liu and Wang2022, Reference Lei, Lu, Gong and Wang2023a
). The generality of this wettability effect transition was demonstrated based on another rock, see Text S6 and figure S6.
Our results suggest that the optimum wetting conditions are closely associated with the geometric confinements of porous media. In 3-D porous media, such as rock and soils, the most favourable wetting conditions for displacement efficiency occur within the weak imbibition regime. Conversely, in 2-D porous media, including microfluidic chips and membranes, the displacement efficiency is higher under strong imbibition conditions. The significant influence of 3-D pore geometries on invasion morphology and subsequently on the
${S_i}{-}{\theta }$
relationship has not been sufficiently addressed in the existing theories of immiscible displacement. Overlooking the complex and synergistic effect of 3-D pore structure and wettability on displacement regimes may obscure the true physics and underlying mechanisms governing multiphase flow dynamics and displacement patterns.
3.2. Pore-scale displacement mechanisms by 3-D pore geometry
To gain comprehensive microscopic mechanisms affected by 3-D pore structure as both wettability and capillary number are altered, we conducted a series of direct numerical simulations over a wide range of
$\theta$
and
$Ca$
in localised regions of both 2-D and 2.5-D porous media. Figure 3(a) illustrates the morphology of multiphase displacement patterns obtained from these simulations at different values of
$\theta$
and
$Ca$
in 2-D porous media. The resulting
$Ca{-}{\theta }$
phase diagram is relatively simple, classifying the patterns into two regions: capillary fingering and compact displacement. In capillary fingering, the morphology reveals a distinct preferential flow pathway, while in compact displacement, a uniform displacement front emerges, characterised by cooperative pore filling across different flow branches. As
$\theta$
decreases, there is a transition from capillary fingering to compact displacement, leading to a more uniform and connected invasion pattern. This monotonic trend is consistent across all values of
$Ca$
, though the critical
$\theta$
for this transition decreases as
$Ca$
increases (figure 3
a).
To provide quantitative insight into these invasion patterns, we calculated two key metrics for each simulation: the saturation of the invading fluid
$S_i$
(figure 3
b) and the Euler number
$\chi$
of the invading phase (figure 3
c) at the breakthrough stage, where
$\chi =\beta _0-\beta _1+\beta _2$
(
$\beta _0$
as the number of connected components,
$\beta _1$
representing the number of redundant loops, and
$\beta _2$
denoting the cavities). A stable and complete displacement typically results in higher
$S_i$
but a smaller
$\chi$
. Negative
$\chi$
indicates connected network structures, signifying complete displacement, while positive values reflect isolated, point-like structures, indicating incomplete displacement. The quantification of
$S_i$
and
$\chi$
follows the monotonic trend, as shown in figures 3(b) and 3(c). For 2-D porous media, we select the Euler number between 30
$^\circ$
and 60
$^\circ$
under
$Ca = 2.6\times {10^{-5}}$
, identified from the phase diagram (figure 3
a), as the critical topological parameter distinguishing stable and unstable displacement (figure 3
c), which is mainly dominated by capillarity for low
$Ca$
.
Invasion phase diagram and quantification results for 2-D porous media under different
$Ca$
by direct numerical simulations. (a) Displacement patterns for different wettability conditions and capillary number. (b) Invading fluid saturation in porous media and (c) Euler number of the invading fluid are functions of the static contact angle
$\theta$
.
However, the typical
${Ca}{-}\theta$
phase diagram is significantly influenced by the 3-D pore geometry, where the displacement regimes can be categorised into four distinct regimes: snap-off-induced trapping, capillary fingering, compact displacement, and coupling of capillary fingering and snap-off (figure 4
a). Under strong imbibition, snap-off occurs due to the instability of the water film flowing along corners of solid surfaces, resulting in the formation of more isolated ganglia. Under strong drainage, the coupling of capillary fingering and snap-off leads to more heterogeneous and discrete flow patterns. Notably, this complex
${Ca}{-}\theta$
phase diagram, with its distinct displacement regimes, was observed for the first time, suggesting that the physical mechanisms governing the two-phase flow in 2-D porous media may not accurately capture all behaviour in 3-D porous media. For 2.5-D porous media, we select the characteristic Euler number between 60
$^\circ$
and 90
$^\circ$
under high
$Ca$
,
$Ca = 2.6\times {10^{-3}}$
, identified from the phase diagram (figure 4
a), as the critical topological parameter distinguishing stable and unstable displacement (figure 4
c), which is mainly controlled by viscous force for high
$Ca$
. The further division of unstable areas is based on qualitative observation.
Invasion phase diagram and quantification results for 2.5-D porous media under different
$Ca$
by direct numerical simulations. (a) Displacement patterns for different wettability conditions and capillary number. Invading fluid saturation in porous media (b) and Euler number of the invading fluid (c) are functions of the static contact angle
$\theta$
.
In comparing the
${Ca}{-}\theta$
phase diagram of 2-D and 2.5-D porous media (figures 3 and 4), the compact displacement at low
$Ca$
and low
$\theta$
in 2-D porous media is replaced by snap-off events in 2.5-D porous media. Similarly, the capillary fingering regime at high
$\theta$
or high
$Ca$
becomes more complex, as the combination of capillary fingering and snap-off in 2.5-D porous media. The critical capillary number for the transition from capillary fingering to the combination of capillary fingering and snap-off is between
$2.6 \times 10^{-4}$
and
$2.6 \times 10^{-3}$
, which is consistent with the critical capillary number proposed by Datta, Dupin & Weitz (Reference Datta, Dupin and Weitz2014) regarding the transition from fully connected to the broken-up flow in 3-D porous media. This comparison provides a valuable reference for using microfluidic experiments to study multiphase displacement mechanisms in 3-D porous media, highlighting the limitations and applicability of 2-D porous media in explaining displacement mechanisms in 3-D porous structures.
Notably, only the transition from snap-off-induced trapping to compact displacement significantly impacts the wettability effect. During drainage, pore-scale events consistently result in lower invading fluid saturation but do not affect the overall wettability trend. A clear non-monotonic wettability rule is observed at low
$Ca$
that as the wettability shifts from strong drainage to strong imbibition, pore-scale events transition from a combination of capillary fingering and snap-off, to pure capillary fingering, and finally to snap-off events at low
$Ca$
. Increasing
$Ca$
suppresses snap-off events, resulting in a monotonic wettability rule, as demonstrated quantitatively by the invading fluid saturation and Euler number (figures 4
b and 4
c). The critical capillary number for the wettability effect transition is approximately
$2.6\times 10^{-4}$
.
3.3. Critical capillary number for non-monotonic wettability effect
To quantify the critical capillary number
$Ca^*$
for the transition from monotonic to non-monotonic wettability effect in porous media with 3-D pore geometry under
$M\gt 1$
, we analyse the balance between viscous and capillary forces during the transition from fully connected to fragmented flow, focusing on the low
$Ca$
regime with small
$\theta$
and the high
$Ca$
regime with large
$\theta$
. The viscous pressure difference along the porous media can be written as (Hilfer et al. Reference Hilfer, Armstrong, Berg, Georgiadis and Ott2015)
${P_v} = {{{\mu _i}\phi {v_i}L} \mathord / k}$
, where
$v_i$
is the velocity of the invading fluid,
$k$
is the permeability of the porous medium,
$\phi$
is the porosity of the porous media, and
$L$
is the length of the porous media. The characteristic capillary pressure in the porous medium is quantified by the Young–Laplace law as
${P_c} = {{2\gamma \cos \theta } \mathord / {{R_c}}}$
, where
$R_c$
is the characteristic radius of flow channels. The transition from stable to unstable flow indicates that
${P_v} \approx {P_c}$
. According to the definition of capillary number, we can get the critical capillary number from
$C{a^*} = {{{\mu _i}{v_i}} \mathord / \gamma } \approx {{2k\cos \theta } \mathord / { ( {\phi L{R_c}} )}}$
, where
$\phi = 0.38$
is the porosity of porous media,
$L=2\,\rm mm$
is the length of porous media,
$R_c=58$
$\unicode{x03BC}$
m is the characteristic hydraulic radius for 2.5-D porous media, the permeability of 2.5-D porous media was calculated by the single-phase network model, and the corresponding permeabilities are
$k=7.1 \times 10^{-12}$
m
$^2$
. Therefore, we can get the corresponding critical capillary number as
$Ca^*\approx 3.8\times 10^{-4}$
, which is qualitatively consistent with the simulation results. Both direct numerical simulations and theoretical predictions highlight that increasing
$Ca$
will lead to the transition from snap-off events to compact displacement and from capillary fingering to a combination of snap-off and capillary fingering. These complex phenomena, triggered by the 3-D pore geometry, dominate the wettability effect during multiphase displacement in porous media. Notably,
$Ca^*$
marks the transition from main meniscus flow to corner or film flow, causing snap-off and incomplete displacement. This differs from the previous critical
$Ca$
, which separates viscous and capillary fingering under
$M \ll 1$
(Holtzman & Segre Reference Holtzman and Segre2015).
4. Conclusion
Through microfluidic experiments, direct simulations and network modelling, we systematically study the transitions of wettability effects during multiphase displacement affected by geometric confinements and capillary numbers. The network modelling and the microfluidic experiments indicate that transitions from monotonic to non-monotonic wettability effects during multiphase displacement can be achieved by manipulating the geometric confinement from 2-D to 3-D porous structures. Microscopic flow dynamics and macroscopic displacement patterns under a wide range of contact angles and capillary numbers were further explored and quantified by direct numerical simulations. For 2-D porous media, a monotonic wettability rule is observed by varying wettability from strong drainage to strong imbibition leading to higher displacement efficiency, where dominated pore-scale events will be transformed from unstable capillary fingering to stable compact displacement. However, for 2.5-D and 3-D porous media, the wettability effect transformed into a non-monotonic rule, which originates from the coupling of capillary fingering and snap-off events under strong drainage and frequent snap-off events under strong imbibition. The 3-D pore geometry triggering frequent snap-off events under strong imbibition leads to the non-monotonic wettability rule. This transition is observed predominantly at lower capillary numbers but is invalid at higher capillary numbers.
However, our results are limited by the assumption of a homogeneous wettability distribution without considering heterogeneous wettability. Additionally, the presented rules apply to a viscously stable case (
$M\gt 1$
), while viscous instability (
$M\ll 1$
) could lead to non-monotonic wettability behaviour (Zhao et al. Reference Zhao, MacMinn and Juanes2016). These issues should be addressed in future research. Nevertheless, we emphasise the limitations and applicability of 2-D porous media models in explaining the displacement mechanisms observed in 3-D porous media. The understanding of wettability effects during multiphase displacement is crucial for controlling multiphase behaviour in naturally occurring porous media, such as soil and rock. Additionally, it provides design principles for artificial porous media (e.g. microfluidic chips, chemical reactors and fuel cell systems) to achieve the desired immiscible displacement behaviour.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2025.10196.
Funding
This work is financially supported by the NSF grant of China (no. 12432013, U24B6003) and the Tsinghua University Initiative Scientific Research Program.
Declaration of interests
The authors report no conflict of interest.
Author contributions
M.W. designed research and supervised the project; W.L. performed research; W.G., X. L., Y. Y. and H.Y. analysed data; and W. L. and M. W. co-wrote the manuscript with feedback from all authors.
