Understanding the solutal convection is a crucial step towards accurate prediction of CO$_2$ sequestration in deep saline aquifers. In this study, pore-scale resolved direct numerical simulations (DNS) are performed to analyse the scaling laws of the solutal convection in porous media. The porous media studied are composed of uniformly distributed square or circular elements. The Rayleigh numbers in the range $426 \le Ra \le 80\,000$, the Darcy numbers in the range $1.7\times 10^{-8} \le Da \le 8.8\times 10^{-6}$ and the Schmidt numbers in the range $250 \le Sc \le 10^4$ are considered in the DNS. The main time, length and velocity scales affecting the solutal convection are classified as the diffusive scales ($t_I$, $l_I$ and $u_I$), the convective scales ($t_{II}$, $l_{II}$ and $u_{II}$) and the shut-down scales ($t_{III}$, $l_{III}$ and $u_{III}$). These scales determine the pore-scale Rayleigh number $Ra_K$ and the Rayleigh number $Ra$. Based on the DNS results, the scaling laws for the transient dissolution flux are proposed in the different regimes of convection. It is found that the onset time of convection ($t_{oc}$) and the period of the flux-growth regime ($\Delta t_{fg}$) vary linearly with the convective time scale $t_{II}$. The merging of the original plumes and the re-initiation of the new plumes start in the same period, resulting in the merging re-initiation regime. Horizontal dispersion plays an important role in plume merging. The dissolution flux $F$ and the time since the onset of convection $t^{\ast }$ have an $F / u_{II} \sim (t^{\ast }/ t_{II})^{-0.2}$ scaling in the later stage of the merging re-initiation regime. This scaling is caused by the merging of the low-wavenumber plumes. It becomes more pronounced with decreasing porosity and leads to the nonlinear relationship between the Sherwood number $Sh$ and $Ra$ when the domain is not high enough for the plumes to fully develop. According to the DNS results, a regime is expected that is independent of both $Ra$ and $Ra_K$, while the dimensionless constant flux $F_{cf}/u_{II}$ in this regime decreases with decreasing porosity. When the mean solute concentration reaches approximately 30 %, convection enters the shut-down regime. For large $Ra$, the dimensionless flux in the shut-down regime follows the scaling law $F/u_{III}\sim (t/t_{III})^{-1.2}$ at large porosity ($\phi =0.56$) and $F/u_{III}\sim (t/t_{III})^{-1.5}$ at small porosity ($\phi =0.36$ or $0.32$). The study shows the significant pore-scale effect on the convection. The DNS cases in this study are mainly for simplified geometries and large $Ra_K$. This can lead to uncertainties of the obtained scaling coefficients. It is important to determine the scaling coefficients for real geological formations to predict a real CO$_2$ sequestration process.

The effect of microbial activity on buoyancy-driven flow within a porous layer is analysed. The input fluid provides an energy source for the growth of biofilms on the porous rock. At each location within the porous layer, the porosity and permeability begin to decrease once the input fluid has invaded. This leads to an evolving rock heterogeneity that depends on the passing time of the input fluid. Hence, the evolution of the flow is partly controlled by its own history. We present an axisymmetric gravity current model, accounting for this effect. In general, a reduction in permeability leads to the flow having a lesser extent in the radial direction and greater thickness (extent in the cross-flow direction), whilst a reduction in porosity has negligible effect on the thickness but leads to a much greater radial extent. The flow is fastest near the free surface where the permeability is greatest. In the case where the porosity and permeability reduce as power-law functions of fluid residence time, the evolution of the flow and the rock properties are self-similar. Consumption of the input fluid by the microbes is also incorporated in the model and it generally leads to flows with lesser radial extent but little change in the thickness. The three impacts of microbial growth (volume loss owing to consumption and the reduction in permeability and porosity) each influence the flow in substantially different ways and the interplay is analysed. A motivation of the study, the underground storage of hydrogen, is briefly discussed.

The present study aims to examine the temporal linear stability analysis of isothermal plane Couette flow over a porous layer using the two-domain approach. The flow in the porous layer is described by the unsteady Darcy–Brinkman equations, whereas it is characterised by the Navier–Stokes equations in the fluid layer. In contrast to the Darcy model, it is observed that the isothermal plane Couette flow becomes unstable for such a superposed system on the inclusion of the Brinkman term. From the stability analysis, the two-dimensional mode is found to be least stable, and two modes of instability, namely porous mode and mixed mode are obtained under the consideration of the Darcy–Brinkman model along with advection term (DBA model). For Darcy number $(\delta )=0.01$, depending on the value of the stress-jump coefficient, mixed mode controls the instability of the system at small values of depth ratio $(\hat {d})$, and it disappears for relatively high values of $\hat {d}$, where the porous mode dominates. In addition, it has been observed that when $\hat {d}=0.1$, the critical mode of instability is found to be mixed for $\delta >0.02$ and porous for $\delta \le 0.02$. The stress-jump coefficient destabilises the flow in terms of energy production through perturbed stresses at the interface. As observed in the case of isothermal plane Poiseuille flow studied by Chang, Chen & Straughan (J. Fluid Mech., vol. 564, 2006, pp. 287–303), here also depth ratio (Darcy number) stabilises (destabilises) the flow. However, this characteristic does not remain valid when the advection term is eliminated from the considered momentum equation. For a certain range of $\hat {d} (\delta )$, the destabilising (stabilising) characteristic of the respective parameters are encountered when the fluid mode of instability prevails.

High-pressure fluid transport in nanoporous media such as shale formations requires further understanding because conventional continuum approaches become inadequate due to their ultralow permeability and complexity of transport mechanisms. We propose a species-based approach for modelling two partially miscible, multicomponent fluids in nanoporous media – one that does not rely on conventional bulk fluid transport frameworks but on species movement. We develop a numerical model for species transport of partially miscible, non-ideal fluid mixtures using the chemical potential gradient as the driving force. The model considers the binary friction concept to include the friction between fluid molecules as well as between fluid molecules and pore walls, and incorporates the key multicomponent transport mechanisms – Knudsen, viscous and molecular diffusion. Under single-phase conditions, the system under consideration is quantified by introducing multicomponent Sherwood number (Sh), Péclet number (Pe) and fluid–solid friction modulus (φ). Despite the complexity of fluid transport in nanopores, the steady-state single-phase transport results reveal the contribution of diffusion in nanopores, where all parameters collapse on a set of master curves for the multicomponent Sh with a dependence on multicomponent Pe and φ. Unsteady state, two-phase transport modelling of the codiffusion process shows that light and intermediate alkanes are produced much higher than heavy alkanes when the vapour phase appears. We demonstrate that the pressure gradient is also crucial in promoting CO2 and alkane mixing during counterdiffusion processes. These results stress the need for a paradigm shift from classical bulk flow modelling to species-based transport modelling in nanoporous media.

We investigate effect of porous insert located upstream of the separation edge of a backward-facing step (BFS) in early transitional regime as a function of Reynolds number. This is an example of hydrodynamic system that is a combination of separated shear flow with large amplification potential and porous materials known for efficient flow destabilisation. Spectral analysis reveals that dynamics of BFS is dominated by spectral modes that remain globally coherent along the streamwise direction. We detect two branches of characteristic frequencies in the flow and with Hilbert transform we characterise their spatial support. For low Reynolds numbers, the dynamics of the flow is dominated by lower frequency, whereas for sufficiently large Reynolds numbers cross-over to higher frequencies is observed. Increasing permeability of the porous insert results in decrease in Reynolds number value, at which frequency cross-over occurs. By comparing normalised frequencies on each branch with local stability analysis, we attribute Kelvin–Helmholtz and Tollmien–Schlichting instabilities to upper and lower frequency branches, respectively. Finally, our results show that porous inserts enhance Kelvin–Helmholtz instability and promote transition to oscillator-type dynamics. Specifically, the amplitude of vortical (BFS) structures associated with higher-frequency branch follows Landau model prediction for all investigated porous inserts.

A macroscopic model for perfect-slip flow in porous media is derived in this work, starting from the pore-scale flow problem and making use of an upscaling technique based on the adjoint method and Green's formula. It is shown that the averaged momentum equation has a Darcy form in which the permeability tensor, $\boldsymbol{\mathsf{K}}_{ps}$, is obtained from an associated adjoint (closure) problem that is to be solved on a (periodic) unit cell representative of the structure. Similarly to the classical permeability tensor, $\boldsymbol{\mathsf{K}}$, in the no-slip regime, $\boldsymbol{\mathsf{K}}_{ps}$ is intrinsic to the porous medium under consideration and is shown to be symmetric and positive. Integral relationships between $\boldsymbol{\mathsf{K}}_{ps}$, the partial-slip flow permeability tensor, $\boldsymbol{\mathsf{K}}_{s}$, and $\boldsymbol{\mathsf{K}}$ are derived. Numerical simulations carried out on two-dimensional model porous structures, together with an approximate analytical solution and an empirical correlation for a particular configuration, confirm the validity of the macroscopic model and the relationship between $\boldsymbol{\mathsf{K}}_{ps}$ and $\boldsymbol{\mathsf{K}}$.

In studying the transport of inclusions in multiphase systems we are often interested in integrated quantities such as the net force and the net velocity of the inclusions. In the reciprocal theorem the known solution to the first and typically easier boundary value problem is used to compute the integrated quantities, such as the net force, in the second problem without the need to solve that problem. Here, we derive a reciprocal theorem for poro-viscoelastic (or biphasic) materials that are composed of a linear compressible solid phase, permeated by a viscous fluid. As an example, we analytically calculate the time-dependent net force on a rigid sphere in response to point forces applied to the elastic network and the Newtonian fluid phases of the biphasic material. We show that when the point force is applied to the fluid phase, the net force on the sphere evolves over time scales that are independent of the distance between the point force and the sphere; in comparison, when the point force is applied to the elastic phase, the time scale for force development increases quadratically with the distance, in line with the scaling of poroelastic relaxation time. Finally, we formulate and discuss how the reciprocal theorem can be applied to other areas, including (i) calculating the network slip on the sphere's surface, (ii) computing the leading-order effects of nonlinearities in the fluid and network forces and stresses, and (iii) calculating self-propulsion in biphasic systems.

The injection of ${\rm CO}_2$ into depleted reservoirs carries the potential for significant Joule–Thomson cooling, when dense, supercritical ${\rm CO}_2$ is injected into a strongly under-pressured reservoir. The resulting low temperatures around the wellbore risk causing thermal fracturing of the well/near-well region or causing freezing of pore waters or formation of gas hydrates which would reduce injectivity and jeopardise well and reservoir integrity. These risks are particularly acute during injection start-up when ${\rm CO}_2$ is in the gas stability field. In this paper we present a model of non-isothermal single-phase flow in the near-wellbore region. We show that during radial injection, with fixed mass injection rate, transient Joule–Thomson cooling can be described by similarity solutions at early times. The positions of the ${\rm CO}_2$ and thermal fronts are described by self-similar scaling relations. We show that, in contrast to steady-state flow, transient flow causes slight heating of ${\rm CO}_2$ and reservoir gas either side of the thermal front, as pressure diffuses into the reservoir. The scaling analysis here identifies the parametric dependence of Joule–Thomson cooling. We present a sensitivity analysis which demonstrates that the primary controls on the degree of cooling are reservoir permeability, reservoir thickness, injection rate and Joule–Thomson coefficient. The analysis presented provides a computationally efficient approach for assessing the degree of Joule–Thomson cooling expected during injection start-up, providing a complement to complex, fully resolved numerical simulations.

We explored the instability dynamics of the viscous fingering interaction in dual displacement fronts by varying the viscosity configuration. Four regimes of rear-dominated fingering, front-dominated fingering, dual fingering and stable were identified. By using the breakthrough time, which refers to the breakup of the dual displacement fronts, the instability dynamics were modelled, and a regime map was developed. These serve as a tool for effectively harnessing the dual displacement fronts for various applications, such as hydrogeology, petroleum, chemical processes and microfluidics.

The nature and behaviour of the drag coefficient $C_D$ of irregularly shaped grains within a wide range of Reynolds numbers $Re$ is discussed. The morphology of the grains is controlled by their fractal description, and they differ in shape. Using computational fluid dynamics tools, the characteristics of the boundary layer at high $Re$ has been determined by applying the Reynolds-averaged Navier–Stokes turbulence model. Both grid resolution and mesh size dependence are validated with well-reported previous experimental results applied in flow around isolated smooth spheres. The drag coefficient for irregularly shaped grains is shown to be higher than that for spherical shapes, also showing a strong drop in its value at high $Re$. This drag crisis is reported at lower $Re$ compared to the smooth sphere, but higher critical $C_D$, demonstrating that the morphology of the particle accelerates this crisis. Furthermore, the dependence of $C_D$ on $Re$ in this type of geometry can be represented qualitatively by four defined zones: subcritical, critical, supercritical and transcritical. The orientational dependence for both particles with respect to the fluid flow is analysed, where our findings show an interesting oscillatory behaviour of $C_D$ as a function of the angle of incidence, fitting the results to a sine-squared interpolation, predicted for particles within the Stokes laminar regime ($Re\ll 1$) and for elongated/flattened spheroids up to $Re=2000$. A statistical analysis shows that this system satisfies a Weibullian behaviour of the drag coefficient when random azimuthal and polar rotation angles are considered.

Previous work suggests that the arrangement of elements in an obstruction may influence the bulk flow velocity through the obstruction, but the physical mechanisms for this influence are not yet clear. This is the motivation for this study, where direct numerical simulation is used to investigate flow through an array of cylinders at a resolution sufficient to observe interactions between wakes of individual elements. The arrangement is altered by varying the gap ratio $G/d$ (1.2 – 18, G is the distance between two adjacent cylinders, d is the cylinder diameter), array-to-element diameter ratio $D/d$ (3.6 – 200, D is the array diameter), and incident flow angle ($0^{\circ} - 30^{\circ}$). Depending on the element arrangement, it is found that the average root-mean-square lift and drag coefficients can vary by an order of magnitude, whilst the average time-mean drag coefficient of individual cylinders ($\overline{C_{d}}$), and the bulk velocity are found to vary by up to $50\,\%$ and a factor of 2, respectively. These arrangement effects are a consequence of the variation in flow and drag characteristics of individual cylinders within the array. The arrangement effects become most critical in the intermediate range of flow blockage parameter $\mathit{\Gamma_{D}^{\prime}} = 0.5-1.5$ ($\mathit{\Gamma_{D}^{\prime}}=\overline{C_{d}}aD/(1-\phi)$, where a is frontal element area per unit volume, and $\phi$ is solid volume fraction), due to the high variability in element-scale flow characteristics. Across the full range of arrangements modelled, it is confirmed that the bulk velocity is governed by flow blockage parameter but only if the drag coefficient incorporates arrangement effects. Using these results, this paper proposes a framework for describing and predicting flow through an array across a variety of arrangements.

Mixing describes the process by which solutes evolve from an initial heterogeneous state to uniformity under the stirring action of a fluid flow. Fluid stretching forms thin scalar lamellae that coalesce due to molecular diffusion. Owing to the linearity of the advection–diffusion equation, coalescence can be envisioned as an aggregation process. Here, we demonstrate that in smooth two-dimensional chaotic flows, mixing obeys a correlated aggregation process, where the spatial distribution of the number of lamellae in aggregates is highly correlated with their elongation, and is set by the fractal properties of the advected material lines. We show that the presence of correlations makes mixing less efficient than a completely random aggregation process because lamellae with similar elongations and scalar levels tend to remain isolated from each other. We show that correlated aggregation is uniquely determined by a single exponent that quantifies the effective number of random aggregation events. These findings expand aggregation theories to a larger class of systems, which have relevance to various fundamental and applied mixing problems.

It is known that the dispersion of colloidal particles in porous media is determined by medium structure, pore-scale flow variability and diffusion. However, much less is known about how diffusiophoresis, that is, the motion of colloidal particles along salt gradients, impacts large-scale particle dispersion in porous media. To shed light on this question, we perform detailed pore-scale simulations of fluid flow, solute transport and diffusiophoretic particle transport in a two-dimensional hyper-uniform porous medium. Particles and solute are initially uniformly distributed throughout the medium. The medium is flushed at constant flow rate, and particle breakthrough curves are recorded at the outlet to assess the macroscopic effects of diffusiophoresis. Particle breakthrough curves show non-Fickian behaviour manifested by strong tailing that is controlled by the diffusiophoretic mobility. Although diffusiophoresis is a short-time, microscopic phenomenon owing to the fast attenuation of salt gradients, it governs macroscopic colloid dispersion through the partitioning of particles into transmitting and dead-end pores. We quantify these behaviours by an upscaled analytical model that describes both the retention and release of colloids in dead-end pores and the observed long-time tailings. Our results suggest that diffusiophoresis is an efficient tool to control particle dispersion and filtration through porous media.

We explore predictions of two models of one-dimensional capillary rise in rigid and partially saturated porous media. One is an existing one from the literature and the second is a free-boundary model based on Richards’ equation with two moving boundaries of the evolving partially saturated region. Both models involve the specification of saturation-dependent functions for local capillary pressure and permeability and connect to classical models for saturated porous media. Existing capillary-rise experiments show two notable regimes: (i) an early-time regime typically well-described by classical capillary-rise theory in a fully saturated porous media, and (ii) a long-time regime that has anomalous dynamics in which the capillary-rise height may scale with a non-classical power law in time or have more complicated dynamics. We demonstrate that the predictions of both models compare well with experimental capillary-rise data over early- and long-time regimes gathered from three independent studies in the literature. The model predictions also shed light on recent scaling laws that relate the capillary pressure and permeability of the partially saturated media to the capillary-rise height. We use these models to probe computationally observed permeability relationships to capillary-rise height. We demonstrate that a recently proposed permeability scaling for the anomalous capillary-rise regime is indeed realized and is particularly apparent in the lower portion of the partially saturated media. For our free-boundary model we also compute capillary pressure measures and show that these reveal the linear relation between the capillary pressure and capillary-rise height expected for a capillarity–gravity balance in the upper portion of the partially saturated porous media.

Utilizing the discrete element method and the pore network model, we numerically investigate the impact of compaction on the longitudinal dispersion coefficient of porous media. Notably, the dispersion coefficient exhibits a non-monotonic dependence on the degree of compaction, which is distinguished by the presence of three distinct regimes in the variation of dispersion coefficient. The non-monotonic variation of dispersion coefficient is attributed to the disparate effect of compaction on dispersion mechanisms. Specifically, the porous medium tightens with an increasing pressure load, reducing the effect of molecular diffusion that primarily governs at small Péclet numbers. On the other hand, heightened pressure loads enhance the heterogeneity of pore structures, resulting in increased disorder and a higher proportion of stagnant zones within porous media flow. These enhancements further strengthen mechanical dispersion and hold-up dispersion, respectively, both acting at higher Péclet numbers. It is crucial to highlight that hold-up dispersion is induced by the low-velocity regions in porous media flow, which differ fundamentally from zero-velocity regions (such as dead-ends or the interior of permeable grains) as described by the classical theory of dispersion. The competition between weakened molecular diffusion and enhanced hold-up dispersion and mechanical dispersion, together with the shift in the dominance of dispersion mechanisms across various Péclet numbers, results in multiple regimes in the variation of dispersion coefficients. Our study provides unique insights into structural design and modulation of the dispersion coefficient of porous materials.

In this study, we propose a novel cooling scheme that utilizes copper inverse opals (CIOs) for surface enhancement in a single-phase impingement jet cooling system. We perform computational fluid dynamics simulations to evaluate the cooling performance of the CIO jet coolers. Our modelling results indicate that the proposed CIO-coated cooler can significantly reduce the average temperature and improve the temperature uniformity across the entire chip surface. The average Nusselt number of the CIO-coated cooler can reach up to 2.8 times that of the flat surface jet cooler. However, the porous structure of the CIO-coated cooler increases the total pressure drop. To determine designs with high cooling performance and low energy consumption, we investigate two crucial design factors, namely the inlet velocity and the nozzle-to-CIO distance. Our analysis reveals that increasing the inlet velocity further enhances the heat transfer, but at the expense of high pressure drop. On the other hand, a larger nozzle-to-CIO distance results in a lower pressure drop but also reduces the heat transfer coefficient. The effects of nozzle-to-CIO distance are further understood by studying the flow resistance network. Furthermore, we present a reduced-order model that accurately captures the thermofluidic characteristics of the proposed design.

A theoretical analysis is presented of peristaltic pumping down a narrow conduit with permeable walls, motivated by the flushing action of lugworms and other marine organisms in sandy burrows. Flow in the conduit is dealt with using lubrication theory; the leakage into the surrounding medium is taken into account by exploiting slender-body theory to solve the associated Darcy problem. By adopting a model for the local force balance on the pumping surface, we bridge between the limits in which the pump operates with either fixed load or displacement. In the latter limit we characterize peristaltic waves with either fixed form or ones that partially collapse the conduit. We construct pump characteristics (the relation between the mean flux and net pressure drop) when the burrow wall is impermeable and pressures are fixed at each end, and compare the results with existing laboratory experiments performed on lugworms. We then consider how the peristaltic dynamics is changed when the wall is made permeable. Last, we consider pumping along an impermeable burrow into a leaky head shaft. The results reveal that the permeability of the conduit wall or end can greatly impact the direction and strength of the recirculating flow.

We consider the process of convective dissolution in a homogeneous and isotropic porous medium. The flow is unstable due to the presence of a solute that induces a density difference responsible for driving the flow. The mixing dynamics is thus driven by a Rayleigh–Taylor instability at the pore scale. We investigate the flow at the scale of the pores using Hele-Shaw type experiment with bead packs, two-dimensional direct numerical simulations and physical models. Experiments and simulations have been specifically designed to mimic the same flow conditions, namely matching porosities, high Schmidt numbers and linear dependency of fluid density with solute concentration. In addition, the solid obstacles of the medium are impermeable to fluid and solute. We characterise the evolution of the flow via the mixing length, which quantifies the extension of the mixing region and grows linearly in time. The flow structure, analysed via the centreline mean wavelength, is observed to grow in agreement with theoretical predictions. Finally, we analyse the dissolution dynamics of the system, quantified through the mean scalar dissipation, and three mixing regimes are observed. Initially, the evolution is controlled by diffusion, which produces solute mixing across the initial horizontal interface. Then, when the interfacial diffusive layer is sufficiently thick, it becomes unstable, forming finger-like structures and driving the system into a convection-dominated phase. Finally, when the fingers have grown sufficiently to touch the horizontal boundaries of the domain, the mixing reduces dramatically due to the absence of fresh unmixed fluid. With the aid of simple physical models, we explain the physics of the results obtained numerically and experimentally. The solute evolution presents a self-similar behaviour, and it is controlled by different length scales in each stage of the mixing process, namely the length scale of diffusion, the pore size and the domain height.

We investigate the three-dimensional (3-D) flow around and through a porous screen for various porosities at high Reynolds number $Re = {O}(10^4)$. Historically, the study of this problem has been focused on two-dimensional cases and for screens spanning completely or partially a channel. Since many recent problems have involved a porous object in a 3-D free flow, we present a 3-D model initially based on Koo & James (J. Fluid Mech., vol. 60, 1973, pp. 513–538) and Steiros & Hultmark (J. Fluid Mech., vol. 853, 2018 pp. 1–11) for screens of arbitrary shapes. In addition, we include an empirical viscous correction factor accounting for viscous effects in the vicinity of the screen. We characterize experimentally the aerodynamic drag coefficient for a porous square screen composed of fibres, immersed in a laminar air flow with various solidities and different angles of attack. We test various fibre diameters to explore the effect of the space between the pores on the drag force. Using PIV and hot wire probe measurements, we visualize the flow around and through the screen, and in particular measure the proportion of fluid that is deviated around the screen. The predictions from the model for drag coefficient, flow velocities and streamlines are in good agreement with our experimental results. In particular, we show that local viscous effects are important: at the same solidity and with the same air flow, the drag coefficient and the flow deviations strongly depend on the Reynolds number based on the fibre diameter. The model, taking into account 3-D effects and the shape of the porous screen, and including an empirical viscous correction factor that is valid for fibrous screens may have many applications including the prediction of water collection efficiency for fog harvesters.

Understanding and controlling fluid entrapment during forced imbibition in porous media is crucial for many natural and industrial applications. However, the microscale physics and macroscopic consequences of fluid entrapment in these geometric-confined porous media remain poorly understood. Here, we introduce a novel multidepth microfluidic chip, which can mitigate the depth confinement of traditional two-dimensional (2-D) microfluidic chips and mimic the wide pore size distribution as natural-occurring three-dimensional (3-D) porous media. Based on microfluidic experiments and direct numerical simulations, we observe the fluid-entrapment scenarios and elucidate the underlying complex interaction between geometric confinement, capillary number and wettability. Increasing depth variation can promote fluid entrapment, whereas increasing capillary number and contact angle yield the opposite effect, which seemingly contradicts conventional expectations in traditional 2-D microfluidic chips. The fluid-entrapment scenario in depth-variable microfluidic chips stems from microscopic interfacial phenomena, classified as snap-off and bypass events. We provide theoretical analyses of these pore-scale events and validate corresponding phase diagrams numerically. It is shown that increasing depth variation triggers snap-off and bypass events. Conversely, a higher capillary number suppresses snap-off events under strong imbibition, and an increased contact angle inhibits bypass events under imbibition. These macroscopic imbibition patterns in microfluidic porous media can be linked with these pore-scale events by improved dynamic pore-network models. Our findings bridge the understanding of forced imbibition between 2-D and 3-D porous media and provide design principles for newly engineered porous media with respect to their desired imbibition behaviours.