<figure><img src="image_1.png" /></figure>

# Wettability impacts residual trapping of immiscible fluids during cyclic injection

Zhongzheng Wang<fnref n="1" />,†, Jean-Michel Pereira<fnref n="2" />, Emilie Sauret<fnref n="1" /> and Yixiang Gan<fnref n="3" />,<fnref n="4" />,†

(Received 28 July 2022; revised 28 February 2023; accepted 12 March 2023)

Understanding the hysteretic behaviour in fluid-fluid displacement processes in porous media is critical in many engineering applications. In this work, we study the quasi-static immiscible displacement process in two-dimensional porous media during cyclic injections in the context of carbon geosequestration. The role of wettability on the residual trapping of CO₂ is investigated numerically using an extended interface tracking algorithm. Despite that higher CO₂ saturation can be achieved in CO₂-wet porous media after the first CO₂ injection, the majority of CO₂ is found to be unstable and can be mobilised during subsequent water injection processes. An improvement in the residual trapping of CO₂ is observed as the number of injection cycles increases, which is associated with the dispersion of continuous CO₂ ganglia into numerous smaller blobs. Compared with either water-wet or CO₂-wet porous media, it is found that less CO₂ is trapped within the neutral-wet ones at equilibrium state after a sufficient number of injection cycles. The hysteretic behaviour of saturation between water/CO₂ injection cycles is found to follow an exponential decay, which eventually reaches a finite value. This process corresponds to the shift of the mobile region during displacement from typical capillary fingering to a less ramified regime, which ultimately converges towards main flow channels. This work highlights the hysteretic behaviour during cyclic injections, providing insights on the wettability impacts on multiphase flow in porous media, which is of great importance in applications such as carbon geosequestration and geological hydrogen storage.

**Key words:** porous media, Hele-Shaw flows, multiphase flow

† Email addresses for correspondence: zhongzheng.wang@qut.edu.au, yixiang.gan@sydney.edu.au

<footnote n="1">¹School of Mechanical, Medical and Process Engineering, Faculty of Engineering, Queensland University of Technology, QLD 4001, Australia</footnote>

<footnote n="2">²Navier, Ecole des Ponts, Univ Gustave Eiffel, CNRS, Marne-la-Vallée, France</footnote>

<footnote n="3">³School of Civil Engineering, The University of Sydney, NSW 2006, Australia</footnote>

<footnote n="4">⁴Sydney Nano, The University of Sydney, NSW 2006, Australia</footnote>
# 1. Introduction

Multiphase flow in porous media has received much attention due to its relevance in various engineering applications, such as enhanced oil recovery, underground hydrogen storage and carbon sequestration (Szulczewski et al. 2012; Lake et al. 2014; Heinemann et al. 2021). During fluid-fluid displacement processes, both fluid properties and flow conditions have been found to strongly impact the invasion process. This was revealed in a phase diagram of displacement patterns in the seminal work by Lenormand, Touboul & Zarcone (1988), where the invasion morphologies including capillary fingering, viscous fingering and stable displacement were demonstrated to be controlled by the viscosity ratio of the two fluids and capillary number, which represents the relative importance of viscous force to capillary force. At the same time, wettability, i.e. the fluids' affinity to the porous medium, has been shown to have a profound influence on the pattern formation of multiphase flow (Trojer, Szulczewski & Juanes 2015; Zhao, MacMinn & Juanes 2016; Blunt 2017; Primkulov et al. 2021; Lei et al. 2022). When the invading fluid is wetting to the porous media, i.e. imbibition processes associated with a contact angle less than 90°, compact growth of invading fluid is expected due to the favoured pore-scale overlap invasion mechanism (Cieplak & Robbins 1988; Holtzman & Segre 2015). However, in a drainage process, i.e. the invading fluid is non-wetting to the solid structure, thin fingers develop during the injection process, which often leads to more trapping of the defending phase.

One interesting question is what happens when the displaced fluid re-enters the porous media, which is of particular importance in the application of geological CO₂ sequestration, one of the promising technologies for mitigating greenhouse gas emissions. In the context of carbon sequestration, during the migration of a supercritical CO₂ plume that is subjected to groundwater flow, water displaces CO₂ at the trailing edge of the plume (Juanes, MacMinn & Szulczewski 2010). This corresponds to the cyclic displacement of immiscible fluids (first CO₂ as invading fluid, then water). In addition, it has been recently shown numerically and experimentally that the residual trapping of CO₂ can be improved via alternating injection of CO₂ and water (Herring, Andersson & Wildenschild 2016; Edlmann et al. 2019; Ahn et al. 2020; Herring et al. 2021). As the number of injection cycles increases, the distribution of CO₂ shifts from large continuous ganglia towards numerous smaller blobs (Ahn et al. 2020; Herring et al. 2021). This morphological evolution of trapped CO₂ not only results in more stable capillary trapping due to a smaller buoyancy force for each individual ganglion, mitigating the risks of leakage, but also facilitates subsequent CO₂ dissolution into water as a result of a greater area-to-volume ratio. However, most studies focusing on cyclic injections of multiphase flow in porous media involve no more than four drainage-imbibition cycles due to limited experimental/computational resources, where the evolution of non-wetting phase trapping cannot be completely captured. Another important aspect that remains relatively unexplored is the role of wettability during cyclic injection. Since a wide range of wetting conditions of porous media for CO₂ storage has been observed (Iglauer, Pentland & Busch 2015), it is vital to understand how wettability affects multiphase flow during cyclic injections.

In the context of carbon sequestration, by further extending a pore-resolved interface tracking algorithm (Wang et al. 2022a), we probe the impact of wettability on the residual trapping of CO₂ during cyclic injection in the quasi-static regime (at vanishing capillary number), where the fluid-fluid interfacial tension dominates the invasion processes. Through extensive numerical simulations, i.e. 15 drainage-imbibition cycles under different contact angles, we highlight the phenomenological signatures of cyclic
injection with focus on the evolution of CO₂ morphology, saturation, and ganglia size and numbers. Then, emphasis is placed on wettability impacts on the hysteretic behaviour of CO₂ saturation. Discussion on the implications of the observations is provided.

## 2. Methods

Recently, a pore-resolved interface tracking algorithm for simulating fluid-fluid displacement processes is proposed with aims to study the impact of particle shape and pore structure on multiphase flow (Wang et al. 2022a,b). The algorithm allows the explicit determination of the critical capillary pressure associated with different pore-scale invasion mechanisms based on the interfacial tension, wettability and pore geometry, and the meniscus with the smallest capillary resistance advances at each computation step. This method is an extension of the model developed by Cieplak & Robbins (1988, 1990), in which the pore-scale mechanisms including burst, touch and overlap are considered for tracking the fluid-fluid interface motion in the pore space. The original method has been successfully applied to reproduce experimental results (Zhao et al. 2014; Hu et al. 2019). However, the recent algorithm is more general in terms of its applicability to porous media with complex geometrical features such as angular particles or irregular channels (Wang, Pereira & Gan 2021; Wang et al. 2022a), as opposed to perfect circular posts in the original model. This is achieved by tracking the node-to-node meniscus movement (as opposed to pore-to-pore jumps in the original model) based on the discretisation of the entire pore space. In our previous work, the algorithm was validated through simulations in a single junction micro-model with different channel widths, wetting conditions and boundary conditions, where we demonstrated that the fluid interface motion is captured accurately, with experimentally observed pore-scale mechanisms reproduced, such as meniscus pinning, cooperative pore filling and corner trapping. Further, the experimentally observed transition from capillary fingering to compact displacement in disordered porous media is also captured by the algorithm as the contact angle varies (Wang et al. 2022a). Detailed description of the algorithm can be found in previous studies (Wang et al. 2022a) and is omitted here for brevity.

In this work, we further extend the recent algorithm for consideration of cyclic injection in porous media. Figure 1(a) shows the final invasion morphology when the invading phase (CO₂) percolates through the porous medium that is initially saturated with the defending fluid (water) with a contact angle of $\theta = 60^{\circ}$ (in this work, the contact angle is measured within water). The pink-to-yellow colour bar represents the invasion sequence. Thus, the fluid displacement corresponds to the first drainage process where the porous medium is non-wetting to the invading fluid (CO₂). The simulation procedure is the same as previous single-time displacement processes and no modification on the algorithm is needed. For the first imbibition (water injection), the status of phases is swapped, i.e. the invading/defending phase is changed to the defending/invading phase. Further, all water ganglia within the porous medium are labelled as trapped and their menisci are deactivated to prevent motion. During the water injection process, if the newly injected water touches an existing water ganglion, the trapping status of all menisci belonging to this ganglion is checked and updated. This process may re-activate menisci, which could allow water to quickly percolate through the porous medium by connecting several ganglia. In summary, the key modification for considering multiphase flow during cyclic injection is to allow the re-activation of menisci that belong to initially isolated invading phase ganglia. The determination of critical capillary pressures for different pore-scale mechanisms remains
<figure><img src="image_2.png" /><figcaption>Figure 1. Interface tracking algorithm for cyclic injection. (a) Invasion morphology at percolation for the first drainage (CO₂ injection) and first imbibition (water injection) processes with water contact angle θ = 60°. (b) Schematics of the successive movements of a meniscus along the surface with a sharp corner, corresponding to unpin events. The effective contact angle at the sharp corner is greater than the intrinsic one. (c) A meniscus touches a grain, splitting into two menisci, corresponding to a touch event. (d) Overlap event, where the two menisci merge into one meniscus. The black dots in the magnified image denote the computation nodes.</figcaption></figure>

unchanged (Wang *et al. 2022a*). Several menisci movement mechanisms are shown in figure 1(b–d).

Non-overlapping circular particles of different sizes have been commonly used as porous media with heterogeneous pore structures (Holtzman & Segre 2015; Trojer *et al.* 2015). These patterned microfluidic flow cells (Hele-Shaw cells) as analogues of natural porous media are simple enough to lead to universal findings and complex enough to have direct relevance to engineered and natural porous media (Zhao *et al.* 2016; Lan *et al.* 2020; Primkulov *et al.* 2021). Here, we consider 3-by-1 rectangular porous media filled with 1330 circular posts arranged on a triangular lattice. The centre-to-centre spacing and average radius are 2.46 and 1, respectively, and a 20% variation in particle size with a uniform distribution is introduced, which leads to an overall porosity of 0.40, which is close to experimentally explored values from several existing works, e.g. 0.38–0.42 (Zhao *et al.* 2014), 0.45 (Zhao *et al.* 2016; Lei *et al.* 2022), 0.55 (Hu *et al.* 2019; Lan *et al.* 2020) and 0.39 (Singh *et al.* 2017). Since the fluid displacement process in the quasi-static regime is considered, the actual value of interfacial tension will not impact the results and it will only scale the capillary pressure. During each simulation, the invading phase (either water or CO₂) is injected from a point inlet at the bottom-centre, and the simulation stops when the invading phase reaches the point outlet, which is located at the top-centre.
Solid walls are imposed at boundaries. As shown in figure 1(a), these boundary conditions correspond to typical core-scale experiments of displacement processes. Gravity is ignored in the simulation as past studies have indicated that the gravity can be insignificant during core-scale experiments (Ruprecht et al. 2014), unlike in field-scale analysis. We note that the consideration of gravity using the interface tracking algorithm is straightforward by incorporating an extra pressure term as a function of fluid density contrast and height. Each injection cycle includes a CO₂ injection process, followed by a water injection process. The injection process corresponds to the controlled-velocity boundary condition with small injection rate, which is the capillary-dominated flow regime with vanishing capillary number. Another type of boundary condition for fluid displacement in porous media is the controlled pressure. Moura et al. (2020) experimentally studied the effects of different boundary conditions during slow drainage processes, and it is found that the invasion morphologies are very similar under these two types of boundary conditions, as long as the flow condition is the capillary-dominated regime. In total, 15 cycles are conducted for contact angle $\theta \in \{45^\circ, 135^\circ\}$ with an increment of 7.5°, covering the wettability range of porous media from water-wet to CO₂-wet. Extreme water-wet ($\theta < 45^\circ$) or CO₂-wet ($\theta > 135^\circ$) porous media are not considered in this study, since the current algorithm does not capture the phenomenon of corner flow in extreme wetting conditions (Zhao et al. 2016) (see Liu et al. (2022) as a recent work on the effect of corner flow under complete wetting condition $\theta = 0^\circ$). The simulations are repeated on five different realisations of porous media with the same statistical parameters, i.e. porosity and particle size distribution. We note that although the results and analysis in this work should be valid for any fluid pairs as long as the injection rate is slow (capillary-dominated regime), we will refer the fluids as water and CO₂ since one of the major applications of cyclic injection is carbon geosequestration.

# 3. Results and discussion

## 3.1. Phenomenological characteristics of cyclic injection

We first present the qualitative observations during cyclic injection processes, focusing on the evolution of fluid distribution morphology, saturation, and ganglia size and numbers. Figure 2 shows the invasion morphologies for the first five injection cycles with water contact angle of 60°. During the first CO₂ injection, i.e. drainage process, the CO₂ (the non-wetting phase, in orange colour) invades the medium, developing thin fingers with trapping of water, a typical signature of capillary fingering. In the first imbibition, water re-enters the domain, occupying the pore space shown in light blue. Thus, the light-blue area also represents the CO₂ that has been displaced out of the domain during the process. We define the active area (or mobile fraction) for a displacement process as the region that has a change of status in phase occupation (orange in drainage, light blue in imbibition). It can be seen that the active area decreases as the cycle number increases (with sequence of injection indicated by black arrows). Also, starting from the first imbibition, the active area can be disconnected, different from a continuous phase as in single-time injection processes. This is due to the re-activation of isolated ganglion when connected to the inlet, allowing jumps of meniscus advancement events. Further, the results indicate that the morphology of the active area tends to converge to the same path as the cycle number increases, regardless of drainage or imbibition.

To more clearly show the evolution of the active area during cyclic injection, the temporal colour-map of the active area is shown in figure 3(a). The pattern of the active area shifts from typical capillary fingering to a less ramified regime, which eventually
<figure><img src="image_3.png" /><figcaption>Figure 2. Invasion morphologies of first five injection cycles with water contact angle $\theta = 60^{\circ}$. Orange represents the newly invaded $CO_2$, which is equivalent to the water that is mobilised during the displacement process. Light blue represents newly invaded water (mobilised $CO_2$). Arrows indicate the sequence.</figcaption></figure>

converges towards a single main flow channel. This behaviour in the mobile region is also observed in recent experiments using a micro-model in the capillary-dominated regime (Ahn et al. 2020). Note that the colour-map in figure 3(a) concerns the first five injection cycles to avoid colour skewness. The change in active area after five cycles is found to be insignificant (supplementary movie 1 is available at https://doi.org/10.1017/jfm.2023.222).

Figure 3(b) plots the $CO_2$ saturation $S_{CO_2}$ during cyclic injection. Red and blue lines, with water contact angle being 60°, correspond to drainage and imbibition processes, respectively. Clearly, hysteretic behaviour in $S_{CO_2}$ is observed with increase/decrease in $S_{CO_2}$ during drainage/imbibition. One parameter of significant importance in the context of carbon geosequestration is the residual $CO_2$ saturation $S_{r,CO_2}$ after water injection, which represents the amount of $CO_2$ that is stably trapped by capillary force and cannot be mobilised by water. The blue-dotted line in figure 3(b) indicates that an increase in $S_{r,CO_2}$ with number of injection cycle is observed, consistent with existing experimental and numerical observations (Edlmann et al. 2019; Ahn et al. 2020; Herring et al. 2021). This implies that the $CO_2$ within the porous medium becomes more stable as the number of injection cycles increases, which can be explained by the evolution of $CO_2$ ganglia size and number, as plotted in figures 3(c) and 3(d), respectively. During water injection (imbibition process indicated by blue lines), the $CO_2$ is broken into more ganglia, which is associated
<figure><img src="image_4.png" /><figcaption>Figure 3. (a) Evolution of the active area during displacement processes. Only the first five cycles are shown to avoid colour-map skewness, as the change in active area is insignificant after five cycles. (b) CO₂ saturation SₐCO₂ as a function of cycle number. The blue-dotted line with circles marks the residual CO₂ saturation SₗCO₂ after water injection. (c) Number of CO₂ ganglia. (d) Normalised average CO₂ ganglia size. The ganglia size is normalised by the average pore size (red region in the inset). In panels (b–d), red-solid and blue-solid lines correspond to drainage and imbibition processes, respectively.</figcaption></figure>

with a reduction in the average ganglia size and increase in ganglia number. However, a decrease in ganglia number and increase in average ganglia size are observed during drainage, corresponding to the newly injected CO₂ connecting isolated CO₂ ganglia during the invasion process. The overall trend of decreasing ganglia size thus leads to the capillary trapping being more stable, which also facilitates subsequent CO₂ dissolution into water as a result of greater area-to-volume ratio. Note that the dissolution of CO₂ is not accounted for in the algorithm, as the focus in the current study is placed on the fluid displacement processes, i.e. residual saturation trapping as opposed to solubility trapping (Matter & Kelemen 2009). It is found that 15 injection cycles are sufficient for the displacement
<figure><img src="image_5.png" /><figcaption>Figure 4. (a) CO₂ saturation as a function of contact angle after: (i) the first CO₂ injection, $S_{r,CO_2}^{N=0}$ (red); (ii) the first water injection, $S_{r,CO_2}^{N=1}$ (blue); and (iii) 15th water injection, $S_{r,CO_2}^{N=15}$ (black). Error bars represent the standard deviation. (b) Difference in saturation between $S_{r,CO_2}^{N=0}$ (red) and $S_{r,CO_2}^{N=15}$, which quantifies the ultimate CO₂ saturation reduction due to cyclic injection.</figcaption></figure>

process to reach the equilibrium state, i.e. no further changes in the saturation hysteresis ($\Delta S_{CO_2}$) or ganglia distribution (supplementary movie 1).

The phenomenological characteristics during cyclic injection reported above for $\theta = 60^{\circ}$, i.e. the hysteretic behaviour in fluid morphology, saturation and ganglia distribution, are also qualitatively observed with different contact angles. In the following, the impacts of wettability during cyclic injection processes are quantitatively analysed.

## 3.2. Effects of wettability on the hysteretic saturation of CO₂

First, the CO₂ saturation after the first CO₂ injection $S_{r,CO_2}^{N=0}$ (after zero injection cycle) is plotted as red line in figure 4(a) with error bars representing the standard deviation of five simulations. Here, $S_{r,CO_2}^{N=0}$ also represents the sweep efficiency of a single-time displacement process, which has been investigated extensively (Trojer et al. 2015; Jung et al. 2016; Zhao et al. 2016; Singh et al. 2017; Wang et al. 2019). The increase in $S_{r,CO_2}^{N=0}$ as the porous medium becomes more wetting to the invading phase was explained by the favoured pore-scale overlap events (Cieplak & Robbins 1988; Holtzman & Segre 2015). Figure 4(a) shows that as the water contact angle increases, i.e. invading phase being more wetting to the porous medium, $S_{r,CO_2}^{N=0}$ increases, consistent with past observations (Jung et al. 2016; Wang et al. 2019, 2022a). Particularly, as the contact angle increases to $\theta = 135^{\circ}$, i.e. in CO₂-wet porous media, a compact CO₂ invasion with almost 100% CO₂ saturation is observed.

Then, we plot the residual CO₂ saturation after the first water injection process, $S_{r,CO_2}^{N=1}$, marking the completion of the first cycle. As shown by the blue line in figure 4(a), a significant decrease in CO₂ saturation can be observed after the water flooding for all wetting conditions, with more significant decrease in CO₂ for greater $\theta$. After 15 cycles, the residual CO₂ saturation $S_{r,CO_2}^{N=15}$, which represents the maximum CO₂ saturation after sufficient numbers of cyclic injections, is shown as a black line. Compared with the residual saturation after only one cycle (blue line), a consistent improvement in CO₂ saturation can be seen, revealing the favourable effect of cyclic injection for enhancing capillary trapping of CO₂. Further, the results also indicate that at an equilibrium state, either water-wet or CO₂-wet porous media are better at trapping CO₂ compared with
weakly water-wet or neutral ones. However, in the context of enhanced oil recovery, these results imply that the neutral wetting condition is ideal compared to either oil-wet or water-wet ones for achieving lowest residual oil saturation, which is consistent with experimental observation in the literature (AlRatrout, Blunt & Bijeljic 2018). It is found that the boundary effect enhances CO₂ trapping near the solid walls for CO₂-wet porous media and suppresses CO₂ trapping for water-wet ones (figure 9a). However, the impacts from the boundary are insignificant and do not influence the conclusion from the results (see Appendix A).

The overall variation in CO₂ saturation under different wetting conditions is quantified by the difference between $S_{r,CO_2}^{N=0}$ and $S_{r,CO_2}^{N=15}$ in figure 4(b), where an increasing trend can be observed. The insets show the CO₂ (white) and water (blue) distribution after the first CO₂ injection (left) and after 15 injection cycles (right) for a water-wet case ($\theta = 45^\circ$) and a CO₂-wet case ($\theta = 135^\circ$), respectively. Under the water-wet condition, the initial injection of CO₂ corresponds to a drainage process, where the invasion morphology falls into the capillary fingering regime described by the theory of invasion percolation (Wilkinson & Willemsen 1983; Lenormand & Zarcone 1985). The resulting CO₂ ganglia tend to reside in relatively large pore spaces, which can be easily trapped during subsequent water flooding processes, thus becoming difficult to mobilise. However, compact distribution of CO₂ after the first injection can be seen when $\theta = 135^\circ$, which is expected due to more pore-scale cooperative pore-filling events (Cieplak & Robbins 1988; Holtzman & Segre 2015), resulting in a displacement efficiency close to 1. Nevertheless, a significant amount of CO₂ can be mobilised during subsequent water flooding processes, dramatically reducing the residual CO₂ saturation in the porous media, which is associated with the large value of ($S_{r,CO_2}^{N=0} - S_{r,CO_2}^{N=15}$). In summary, although a CO₂-wet porous medium can have larger capacity in storing CO₂ during the first CO₂ injection, a significant amount of CO₂ is unstable and can be mobilised by subsequent underground water flow, increasing the risk of leakage and leading to a reduced storage capacity of the geological formation.

Another common approach to represent the hysteretic behaviour during cyclic injection is through the capillary pressure-saturation curves. The dimensionless capillary pressure is calculated by $P_c^* = 1/r^*$, where $r^*$ is the dimensionless meniscus radius. Due to the relatively uniform grain size distribution in the absence of gravity, the capillary pressure signal is found to be flat (figure 10). For the ease of visualisation, the capillary pressure is converted by constraining a positive increment in $P_c^*$, such that $P_c^*$ monotonically increases with the invading phase saturation (see Appendix B). Figure 5 shows the history of pressure-saturation relation (averaged from five simulations) during the first five injection cycles for contact angles $\theta = \{45^\circ, 90^\circ, 135^\circ\}$. The increase in the cycle number is denoted by the arrow direction as well as the increasing colour intensity. In water-wet porous media ($\theta = 45^\circ$), it can be seen that the decrease in CO₂ saturation hysteresis during the drainage–imbibition cycles is mainly due to the increase in the residual CO₂ saturation after water flooding, implying that a significant amount of CO₂ is broken up by injected water into smaller blobs during the imbibition process, which leads to more stable capillary trapping. For CO₂-wet porous media ($\theta = 135^\circ$), in contrast, the decrease in saturation hysteresis results from the decrease in CO₂ saturation after CO₂ injection. This is associated with an increasing amount of water being stably trapped within the porous media, reducing the pore space available to CO₂. For neutral-wet porous media, the saturation hysteresis decreases during both water injection and CO₂ injection, which is reflected by the relatively symmetrical narrowing of saturation towards the middle. In all plots, it can be seen that the change in $P_c^*$ is not significant during the injection process due to relatively uniform pore size distribution (20% variation in particle size). Also, vertical
<figure><img src="image_6.png" /><figcaption>Figure 5. Capillary pressure versus CO₂ saturation curves for the first five cycles with contact angles $\theta = \{45^{\circ}, 90^{\circ}, 135^{\circ}\}$ (from left to right). Brown and blue colours respectively represent the CO₂ and water injections. Dashed arrows as well as the increasing colour intensity indicate the increasing cycle number.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. (a) Evolution of CO₂ saturation hysteresis $\Delta S_{CO_2}$ with number of injection cycles N for different contact angles. Solid curves are from (3.1). (b) Fitting parameter a in (3.1) as a function of contact angle. Error bars represent standard deviation.</figcaption></figure>

jumps in $P_c^*$ can be observed between each injection process, indicating a sudden shift in the wettability of the invading fluid at the inlet. Note that the pressure scale in the y-axis for neutral-wet porous media is much smaller, which is expected as both water and CO₂ in this case have the same contact angle $\theta = 90^{\circ}$.

To evaluate how the saturation hysteresis $\Delta S_{CO_2}$ evolves during cyclic injection, figure 6(a) plots $\Delta S_{CO_2}$ against the number of injection cycles with the experimental data from Ahn et al. (2020) added for comparison. It is found that the data can be described by the following exponential decay:

$$ \Delta S_{CO_2} = (S_{r,CO_2}^{N=0} - \Delta S_{CO_2}^*) \cdot a^{-N} + \Delta S_{CO_2}^*, \quad (3.1) $$

with N the number of cycles, a a fitting parameter, $S_{r,CO_2}^{N=0}$ the initial CO₂ saturation after the first CO₂ injection and $\Delta S_{CO_2}^* = \Delta S_{CO_2}^{N=15}$ the saturation hysteresis after a sufficient number of cycles (15 in this study) for the system to reach equilibrium, i.e. the terminal/equilibrium saturation hysteresis. The fitted curves are plotted as solid lines in figure 6(a); they can well capture the evolution of $\Delta S_{CO_2}$. It is found that there is not a significant dependence of a, the decay rate of saturation hysteresis, on contact angles (figure 6b). It is likely that the value of a depends on factors such as connectivity of pore space and particle shape, which is worth investigating in a future study.
Wettability impacts residual trapping of immiscible fluids

<figure><img src="image_8.png" /><figcaption>Figure 7. Evolution of initial and residual CO₂ saturation during cyclic injections for different contact angles. Dashed arrows indicate the direction of increasing cycles for θ = {45°, 90°, 135°}. Black-solid curves correspond to the Land model with parameter C = {0.4, 2}.</figcaption></figure>

We then plot the initial and residual CO₂ saturation over the course of 15 cyclic injections in figure 7 to compare with the classic trapping model by Land (1968). Data from two recent experimental works on cyclic injection are also added (Ahn et al. 2020; Herring et al. 2021). The direction of increasing cycles are marked by dashed arrows. The movement in the initial-residual CO₂ saturation map represents the evolution of saturation hysteresis during cyclic injections. It can be seen that the trajectories from both the current work and literature do not strictly follow the Land model, i.e. $S_r = S_i/(1 + CS_i)$ (Land 1968) (black-solid curves). However, one common feature is that they all travel towards a direction with decreasing C, which is associated with smaller changes in saturation after water flooding. This can be graphically represented by the decrease in the distance with the black-dashed line (zero hysteresis). Further, the stabilisation/convergence of the trajectories towards a certain point provides evidence that the system has reached the hysteresis equilibrium, i.e. the variation in CO₂ saturation between injection cycles becomes constant and does not further change as injection cycle increases. This is however not captured in the past experimental data due to a limited number of injection cycles.

### 3.3. Terminal saturation hysteresis

It can be seen in figures 5 and 6(a) that the hysteresis in saturation does not further change after a sufficient number of cycles, i.e. when the quasi-static cyclic injection processes reach the equilibrium. The non-zero residual hysteresis in saturation $\Delta S_{CO_2}^* = \Delta S_{CO_2}^{N=15}$, or the terminal saturation hysteresis, highlights the history dependence of the quasi-static fluid displacement processes, in line with recent experimental work (Holtzman et al. 2020). The terminal saturation hysteresis $\Delta S_{CO_2}^*$ as a function of contact angle is plotted in figure 8. The $\Delta S_{CO_2}^*$ also represents the converged active area, which is shown in the insets of figure 8 for sample cases with $\theta = \{45°, 90°, 135°\}$. It can be seen that $\Delta S_{CO_2}^*$ increases with contact angle, reaching the peak values at neutral-wet porous media which is associated with the formation of several loops in the active area, before dropping as $\theta$ further increases. The interesting symmetry in $\Delta S_{CO_2}^*$ can be explained by the symmetrical
<figure><img src="image_9.png" /><figcaption>Figure 8. Terminal saturation hysteresis $\Delta S^{*}_{CO_2}$ for different contact angles. Error bars represent standard deviation. Insets show examples of the active area during cyclic injections for contact angles $\theta = \{45^{\circ}, 90^{\circ}, 135^{\circ}\}$.</figcaption></figure>

wettability during cyclic injection, where the contact angle is reversed between cycles. The greater saturation hysteresis, or active area as shown in the insets, under neutral-wet condition is consistent with a recent study where a higher relative permeability is found in neutral-wet porous media compared to water-wet ones (Hashemi, Blunt & Hajibeygi 2021).

To interpret the hysteresis during immiscible fluid displacement in porous media, a discrete-domain model was recently proposed by Cueto-Felgueroso & Juanes (2016) where each sub-domain (REV) of the porous media is represented by a capillary tube whose diameter varies non-monotonically with the axial position, such that the pressure-saturation relationship is also non-monotonic. According to the model, the macroscopically observed hysteresis is a result of the collective behaviour of these interconnected multistable capillary tubes. Thus, the magnitude of the hysteresis should depend on the local energy barrier of the capillary tube when transitioning across different local minima of the Gibbs-like free energy. In the case (current work) of porous media of pore structures with statistically fixed geometrical features, i.e. fixed fluctuations in the pore/throat sizes, the fluid-fluid interfaces in either water-wet ($\theta < 90^{\circ}$) or CO$_2$-wet ($\theta > 90^{\circ}$) will be associated with greater absolute capillary pressure compared with neutral-wet media ($\theta \sim 90^{\circ}$), assuming $P^{*}_{c} \propto 1/\cos\theta$ in capillary tubes (Blunt 2017). This implies that the associated energy barrier between different meta-stable states is also magnified, making the menisci more stable and less likely to be mobilised, consequently leading to less saturation hysteresis between injection cycles, in line with our results in figure 8.

In this work, we have investigated the trapping of immiscible fluids during cyclic injections in the quasi-static regime. In the context of carbon geosequestration, our focus has been placed on investigating the residual saturation trapping mechanism (capillary trapping). We want to point out that if other factors such as viscosity, dissolution or wettability alteration are taken into account, different behaviours in the saturation hysteresis may be expected. It has been established that the trapped phase can be mobilised via momentum transfer (Zarikos et al. 2018) and mechanisms such as droplet fragmentation can occur under high viscous stress (Pak et al. 2015). Furthermore, variations in solid surface wetting condition are possible due to long-term exposure to the
non-wetting phase, high temperature or pressure (Chiquet, Broseta & Thibeau 2007; Kim et al. 2012; Wan, Kim & Tokunaga 2014), which may affect the subsequent displacement processes during cyclic injections. In addition, the heterogeneity of the porous media in the current work is mainly from the variations in particle size. Other forms of heterogeneity, such as wettability (Zhao et al. 2018; Jahanbakhsh, Shahrokhi & Maroto-Valer 2021; Irannezhad et al. 2023), grain shape (Rokhforouz & Akhlaghi Amiri 2019; Wang et al. 2021) and surface roughness (Tanino, Ibekwe & Pokrajac 2020; Wang, Pereira & Gan 2020; Zulfiqar et al. 2020), can also impact the immiscible fluid displacement processes and residual trapping. The results from the current work should lay the foundation for future investigations on the influence of these complexities during cyclic injections.

## 4. Conclusions

By extending a recently developed interface tracking algorithm, the cyclic immersible fluid displacement process is investigated. We present, for the first time to the best of our knowledge, the complete evolution of invasion morphology and saturation change up to the displacement process reaching the equilibrium state, which is associated with constant residual saturation and saturation hysteresis. In the context of carbon geosequestration, the results show that the increase in trapped CO₂ saturation after water injection during cyclic injection, i.e. the increased stability of trapped CO₂, is correlated with an increase in ganglia number and decrease in ganglia size.

We then focus on quantitative characterisation of wettability impacts on the hysteretic behaviour of CO₂ saturation during cyclic injection. Despite that higher CO₂ saturation can be achieved in CO₂-wet porous media after the first CO₂ injection, the majority of CO₂ is unstable and can be mobilised during the water injection process. As the number of injection cycles increases, it is shown that the residual CO₂ saturation after water flooding is consistently improved regardless of the wetting condition. Through analysing the evolution of saturation hysteresis with injection cycle, it is found that the saturation hysteresis follows an exponential decay. Compared with either water-wet or CO₂-wet porous media, less CO₂ is trapped within the neutral-wet ones at equilibrium state. In addition, the terminal CO₂ saturation hysteresis, quantified by the saturation change between water/CO₂ injections after sufficient cycles, is greater in neutral-wet porous media than either water-wet or CO₂-wet ones. Our results compare favourably with the existing, albeit limited amount of, experimental data on cyclic injection, and the insights gained in the current work could help deepen the understanding of multiphase flow during cyclic injections, which is of great importance in engineering applications such as carbon geosequestration and geological hydrogen storage.

**Supplementary movie.** Supplementary movie is available at https://doi.org/10.1017/jfm.2023.222.

**Declaration of interest.** The authors report no conflict of interest.

**Funding.** This work was financially supported by The University of Sydney SOAR Fellowship. Y.G. acknowledges the financial support of Labex MMCD(ANR-11-LABX-022-01) for his stay at Laboratoire Navier at ENPC.

**Author ORCIDs.**

© Zhongzheng Wang https://orcid.org/0000-0001-9456-4744;

© Jean-Michel Pereira https://orcid.org/0000-0002-0290-5191;

© Emilie Sauret https://orcid.org/0000-0002-8322-3319;

© Yixiang Gan https://orcid.org/0000-0002-9621-0277.
<figure><img src="image_10.png" /><figcaption>Figure 9. Boundary effects on the CO₂ distribution. (a) CO₂ distribution along the normalised domain width for contact angles θ = {45° 90°, 135°}. The insets show the raw data (periodically fluctuating curves) and the moving average with span length of one pore size. (b) CO₂ saturation at different injection cycles (same notation as figure 4a) calculated from the area between 0.2 and 0.8 normalised domain width.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10. (a) An example of the raw pressure signal for θ = 90° during CO₂ injection. (b) Magnified region of the first 300 computation steps. Red circles mark the sampled pressure with a positive increment.</figcaption></figure>

**Appendix A. Boundary effect**

To evaluate the sidewall (longer side) boundary effect on the CO₂ distribution, the CO₂ distribution after 15 injection cycles along the horizontal direction is calculated by summing all the CO₂ occupied area (pixels) along the vertical direction. The inset in figure 9(a) shows the raw data (periodically fluctuating line) and the moving average (with span length of one pore size) for θ = 45°. The periodic fluctuation is attributed to the regular (triangular) placement of particles. Figure 9(a) reveals that for water-wet to neutral porous media, the CO₂ occupation of pore space is less favourable along the boundary, whereas for CO₂ wet porous media, more CO₂ is distributed along the boundary walls. To exclude the boundary wall effect, the residual CO₂ saturation at different injection cycles is plotted again, but only with domain region between 0.2 and 0.8 normalised width considered (figure 9b). Compared with figure 4(a), the same qualitative results on CO₂ can
be observed, implying the boundary effect should not influence the conclusions derived
from the main text.

**Appendix B. Pressure signal sampling and conversion**

An example of the raw capillary pressure signal with θ = 90° during the first CO₂ injection process is shown in figure 10(a). Due to the relatively uniform pore size distribution, i.e. 20% variation in particle size, the change in P*ₐ is not significant during the injection process. Figure 10(b) shows the sampling method of pressure signals plotted in figure 5: starting with the initial pressure, the subsequent pressure and associated computation step is only recorded if it is greater than the previous one.

REFERENCES

AHN, H., KIM, S.-O., LEE, M. & WANG, S. 2020 Migration and residual trapping of immiscible fluids during cyclic injection: pore-scale observation and quantitative analysis. *Geofluids* **2020**, 4569208.

ALRATROUT, A., BLUNT, M.J. & BIJELJIC, B. 2018 Wettability in complex porous materials, the mixed-wet state, and its relationship to surface roughness. *Proc. Natl Acad. Sci.* **115** (36), 8901–8906.

BLUNT, M.J. 2017 *Multiphase Flow in Permeable Media: A Pore-Scale Perspective*. Cambridge University Press.

CHIQUET, P., BROSETA, D. & THIBEAU, S. 2007 Wettability alteration of caprock minerals by carbon dioxide. *Geofluids* **7** (2), 112–122.

CIEPLAK, M. & ROBBINS, M.O. 1988 Dynamical transition in quasistatic fluid invasion in porous media. *Phys. Rev. Lett.* **60**, 2042–2045.

CIEPLAK, M. & ROBBINS, M.O. 1990 Influence of contact angle on quasistatic fluid invasion of porous media. *Phys. Rev. B* **41**, 11508–11521.

CUETO-FELGUEROSO, L. & JUANES, R. 2016 A discrete-domain description of multiphase flow in porous media: rugged energy landscapes and the origin of hysteresis. *Geophys. Res. Lett.* **43** (4), 1615–1622.

EDLMANN, K., HINCHLIFFE, S., HEINEMANN, N., JOHNSON, G., ENNIS-KING, J. & MCDERMOTT, C.I. 2019 Cyclic CO₂ – H₂O injection and residual trapping: implications for CO₂ injection efficiency and storage security. *Intl J. Greenhouse Gas Control* **80**, 1–9.

HASHEMI, L., BLUNT, M. & HAJIBEYGI, H. 2021 Pore-scale modelling and sensitivity analyses of hydrogen-brine multiphase flow in geological porous media. *Sci. Rep.* **11** (1), 8348.

HEINEMANN, N., et al. 2021 Enabling large-scale hydrogen storage in porous media – the scientific challenges. *Energy Environ. Sci.* **14**, 853–864.

HERRING, A.L., ANDERSSON, L. & WILDENSCHILD, D. 2016 Enhancing residual trapping of supercritical CO₂ via cyclic injections. *Geophys. Res. Lett.* **43** (18), 9677–9685.

HERRING, A.L., SUN, C., ARMSTRONG, R.T., LI, Z., MCCLURE, J.E. & SAADATFAR, M. 2021 Evolution of bentheimer sandstone wettability during cyclic scCO₂-brine injections. *Water Resour. Res.* **57** (11), e2021WR030891.

HOLTZMAN, R., DENTZ, M., PLANET, R. & ORTÍN, J. 2020 The origin of hysteresis and memory of two-phase flow in disordered media. *Commun. Phys.* **3** (1), 222.

HOLTZMAN, R. & SEGRE, E. 2015 Wettability stabilizes fluid invasion into porous media via nonlocal, cooperative pore filling. *Phys. Rev. Lett.* **115**, 164501.

HU, R., LAN, T., WEI, G.-J. & CHEN, Y.-F. 2019 Phase diagram of quasi-static immiscible displacement in disordered porous media. *J. Fluid Mech.* **875**, 448–475.

IGLAUER, S., PENTLAND, C.H. & BUSCH, A. 2015 CO₂ wettability of seal and reservoir rocks and the implications for carbon geo-sequestration. *Water Resour. Res.* **51** (1), 729–774.

IRANNEZHAD, A., PRIMKULOV, B.K., JUANES, R. & ZHAO, B. 2023 Fluid-fluid displacement in mixed-wet porous media. *Phys. Rev. Fluids* **8**, L012301.

JAHANBAKHSH, A., SHAHROKHI, O. & MAROTO-VALER, M.M. 2021 Understanding the role of wettability distribution on pore-filling and displacement patterns in a homogeneous structure via quasi 3D pore-scale modelling. *Sci. Rep.* **11** (1), 17847.

JUANES, R., MACMINN, C.W. & SZULCZEWSKI, M.L. 2010 The footprint of the CO₂ plume during carbon dioxide storage in saline aquifers: storage efficiency for capillary trapping at the basin scale. *Transp. Porous Med.* **82** (1), 19–30.
JUNG, M., BRINKMANN, M., SEEMANN, R., HILLER, T., SANCHEZ DE LA LAMA, M. & HERMINGHAUS, S. 2016 Wettability controls slow immiscible displacement through local interfacial instabilities. *Phys. Rev. Fluids* 1, 074202.

KIM, Y., WAN, J., KNEAFSEY, T.J. & TOKUNAGA, T.K. 2012 Dewetting of silica surfaces upon reactions with supercritical CO₂ and brine: pore-scale studies in micromodels. *Environ. Sci. Technol.* 46 (7), 4228–4235.

LAKE, L.W., JOHNS, R., ROSSEN, B. & POPE, G. 2014 *Fundamentals of Enhanced Oil Recovery*. Society of Petroleum Engineers.

LAN, T., HU, R., YANG, Z., WU, D.-S. & CHEN, Y.-F. 2020 Transitions of fluid invasion patterns in porous media. *Geophys. Res. Lett.* 47 (20), e2020GL089682.

LAND, C.S. 1968 Calculation of imbibition relative permeability for two- and three-phase flow from rock properties. *Soc. Petrol. Engng J.* 8 (02), 149–156.

LEI, W., LU, X., LIU, F. & WANG, M. 2022 Non-monotonic wettability effects on displacement in heterogeneous porous media. *J. Fluid Mech.* 942, R5.

LENORMAND, R., TOUBOUL, E. & ZARCONE, C. 1988 Numerical models and experiments on immiscible displacements in porous media. *J. Fluid Mech.* 189, 165–187.

LENORMAND, R. & ZARCONE, C. 1985 Invasion percolation in an etched network: measurement of a fractal dimension. *Phys. Rev. Lett.* 54, 2226–2229.

LIU, Y., BERG, S., JU, Y., WEI, W., KOU, J. & CAI, J. 2022 Systematic investigation of corner flow impact in forced imbibition. *Water Resour. Res.* 58 (10), e2022WR032402.

MATTER, J.M. & KELEMEN, P.B. 2009 Permanent storage of carbon dioxide in geological reservoirs by mineral carbonation. *Nat. Geosci.* 2 (12), 837–841.

MOURA, M., MÅLØY, K.J., FLEKKØY, E.G. & TOUSSAINT, R. 2020 Intermittent dynamics of slow drainage experiments in porous media: characterization under different boundary conditions. *Front. Phys.* 7, 217.

PAK, T., BUTLER, I.B., GEIGER, S., VAN DIJKE, M.I.J. & SORBIE, K.S. 2015 Droplet fragmentation: 3D imaging of a previously unidentified pore-scale process during multiphase flow in porous media. *Proc. Natl Acad. Sci.* 112 (7), 1947–1952.

PRIMKULOV, B.K., PAHLAVAN, A.A., FU, X., ZHAO, B., MACMINN, C.W. & JUANES, R. 2021 Wettability and Lenormand’s diagram. *J. Fluid Mech.* 923, A34.

ROKHFOROUZ, M.R. & AKHLAGHI AMIRI, H.A. 2019 Effects of grain size and shape distribution on pore-scale numerical simulation of two-phase flow in a heterogeneous porous medium. *Adv. Water Resour.* 124, 84–95.

RUPRECHT, C., PINI, R., FALTA, R., BENSON, S. & MURDOCH, L. 2014 Hysteretic trapping and relative permeability of CO₂ in sandstone at reservoir conditions. *Intl J. Greenhouse Gas Control* 27, 15–27.

SINGH, K., SCHOLL, H., BRINKMANN, M., MICHIEL, M.D., SCHEEL, M., HERMINGHAUS, S. & SEEMANN, R. 2017 The role of local instabilities in fluid invasion into permeable media. *Sci. Rep.* 7 (1), 444.

SZULCZEWSKI, M.L., MACMINN, C.W., HERZOG, H.J. & JUANES, R. 2012 Lifetime of carbon capture and storage as a climate-change mitigation technology. *Proc. Natl Acad. Sci.* 109 (14), 5185–5189.

TANINO, Y., IBEKWE, A. & POKRAJAC, D. 2020 Impact of grain roughness on residual nonwetting phase cluster size distribution in packed columns of uniform spheres. *Phys. Rev. E* 102, 013109.

TROJER, M., SZULCZEWSKI, M.L. & JUANES, R. 2015 Stabilizing fluid-fluid displacements in porous media through wettability alteration. *Phys. Rev. Appl.* 3, 054008.

WAN, J., KIM, Y. & TOKUNAGA, T.K. 2014 Contact angle measurement ambiguity in supercritical CO₂-water-mineral systems: mica as an example. *Intl J. Greenhouse Gas Control* 31, 128–137.

WANG, Z., CHAUHAN, K., PEREIRA, J.-M. & GAN, Y. 2019 Disorder characterization of porous media and its effect on fluid displacement. *Phys. Rev. Fluids* 4, 034305.

WANG, Z., PEREIRA, J.-M. & GAN, Y. 2020 Effect of wetting transition during multiphase displacement in porous media. *Langmuir* 36 (9), 2449–2458.

WANG, Z., PEREIRA, J.-M. & GAN, Y. 2021 Effect of grain shape on quasi-static fluid-fluid displacement in porous media. *Water Resour. Res.* 57 (4), e2020WR029415.

WANG, Z., PEREIRA, J.-M., SAURET, E., ARYANA, S.A., SHI, Z. & GAN, Y. 2022a A pore-resolved interface tracking algorithm for simulating multiphase flow in arbitrarily structured porous media. *Adv. Water Resour.* 162, 104152.

WANG, Z., PEREIRA, J.-M., SAURET, E. & GAN, Y. 2022b Emergence of unstable invasion during imbibition in regular porous media. *J. Fluid Mech.* 941, A40.

WILKINSON, D & WILLEMSEN, J.F. 1983 Invasion percolation: a new form of percolation theory. *J. Phys. A: Math. Gen.* 16 (14), 3365–3376.
ZARIKOS, I., TERZIS, A., HASSANIZADEH, S. & WEIGAND, B. 2018 Velocity distributions in trapped and mobilized non-wetting phase ganglia in porous media. *Sci. Rep.* **8**, 13228.

ZHAO, J., KANG, Q., YAO, J., VISWANATHAN, H., PAWAR, R., ZHANG, L. & SUN, H. 2018 The effect of wettability heterogeneity on relative permeability of two-phase flow in porous media: a lattice Boltzmann study. *Water Resour. Res.* **54** (2), 1295–1311.

ZHAO, B., MACMINN, C.W., HUPPERT, H.E. & JUANES, R. 2014 Capillary pinning and blunting of immiscible gravity currents in porous media. *Water Resour. Res.* **50** (9), 7067–7081.

ZHAO, B., MACMINN, C.W. & JUANES, R. 2016 Wettability control on multiphase flow in patterned microfluidics. *Proc. Natl Acad. Sci.* **113** (37), 10251–10256.

ZULFIQAR, B., VOGEL, H., DING, Y., GOLMOHAMMADI, S., KÜCHLER, M., REUTER, D. & GEISTLINGER, H. 2020 The impact of wettability and surface roughness on fluid displacement and capillary trapping in 2-D and 3-D porous media: 2. combined effect of wettability, surface roughness, and pore space structure on trapping efficiency in sand packs and micromodels. *Water Resour. Res.* **56** (10), e2020WR027965.