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 geosequestration (Szulczewski et al. Reference Szulczewski, MacMinn, Herzog and Juanes2012; Lake et al. Reference Lake, Johns, Rossen and Pope2014; Heinemann et al. Reference Heinemann2021). 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 (Reference Lenormand, Touboul and Zarcone1988), 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 Reference Trojer, Szulczewski and Juanes2015; Zhao, MacMinn & Juanes Reference Zhao, MacMinn and Juanes2016; Blunt Reference Blunt2017; Primkulov et al. Reference Primkulov, Pahlavan, Fu, Zhao, MacMinn and Juanes2021; Lei et al. Reference Lei, Lu, Liu and Wang2022). When the invading fluid is wetting to the porous media, i.e. imbibition processes associated with a contact angle less than 
$90^\circ$, compact growth of invading fluid is expected due to the favoured pore-scale overlap invasion mechanism (Cieplak & Robbins Reference Cieplak and Robbins1988; Holtzman & Segre Reference Holtzman and Segre2015). 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
$_2$ sequestration, one of the promising technologies for mitigating greenhouse gas emissions. In the context of carbon sequestration, during the migration of a supercritical CO
$_2$ plume that is subjected to groundwater flow, water displaces CO
$_2$ at the trailing edge of the plume (Juanes, MacMinn & Szulczewski Reference Juanes, MacMinn and Szulczewski2010). This corresponds to the cyclic displacement of immiscible fluids (first CO
$_2$ as invading fluid, then water). In addition, it has been recently shown numerically and experimentally that the residual trapping of CO
$_2$ can be improved via alternating injection of CO
$_2$ and water (Herring, Andersson & Wildenschild Reference Herring, Andersson and Wildenschild2016; Edlmann et al. Reference Edlmann, Hinchliffe, Heinemann, Johnson, Ennis-King and McDermott2019; Ahn et al. Reference Ahn, Kim, Lee and Wang2020; Herring et al. Reference Herring, Sun, Armstrong, Li, McClure and Saadatfar2021). As the number of injection cycles increases, the distribution of CO
$_2$ shifts from large continuous ganglia towards numerous smaller blobs (Ahn et al. Reference Ahn, Kim, Lee and Wang2020; Herring et al. Reference Herring, Sun, Armstrong, Li, McClure and Saadatfar2021). This morphological evolution of trapped CO
$_2$ 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
$_2$ 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
$_2$ storage has been observed (Iglauer, Pentland & Busch Reference Iglauer, Pentland and Busch2015), 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. Reference Wang, Pereira, Sauret, Aryana, Shi and Gan2022a), we probe the impact of wettability on the residual trapping of CO
$_2$ 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
$_2$ morphology, saturation, and ganglia size and numbers. Then, emphasis is placed on wettability impacts on the hysteretic behaviour of CO
$_2$ 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. Reference Wang, Pereira, Sauret, Aryana, Shi and Gan2022a,Reference Wang, Pereira, Sauret and Ganb). 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 (Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990), 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. Reference Zhao, MacMinn, Huppert and Juanes2014; Hu et al. Reference Hu, Lan, Wei and Chen2019). 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 Reference Wang, Pereira and Gan2021; Wang et al. Reference Wang, Pereira, Sauret, Aryana, Shi and Gan2022a), 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. Reference Wang, Pereira, Sauret, Aryana, Shi and Gan2022a). Detailed description of the algorithm can be found in previous studies (Wang et al. Reference Wang, Pereira, Sauret, Aryana, Shi and Gan2022a) 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
$_2$) 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
$_2$). 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 unchanged (Wang et al. Reference Wang, Pereira, Sauret, Aryana, Shi and Gan2022a). Several menisci movement mechanisms are shown in figure 1(b–d).
Interface tracking algorithm for cyclic injection. (a) Invasion morphology at percolation for the first drainage (CO
$_2$ injection) and first imbibition (water injection) processes with water contact angle 
$\theta =60^\circ$.(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.
Non-overlapping circular particles of different sizes have been commonly used as porous media with heterogeneous pore structures (Holtzman & Segre Reference Holtzman and Segre2015; Trojer et al. Reference Trojer, Szulczewski and Juanes2015). 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. Reference Zhao, MacMinn and Juanes2016; Lan et al. Reference Lan, Hu, Yang, Wu and Chen2020; Primkulov et al. Reference Primkulov, Pahlavan, Fu, Zhao, MacMinn and Juanes2021). 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. Reference Zhao, MacMinn, Huppert and Juanes2014), 0.45 (Zhao et al. Reference Zhao, MacMinn and Juanes2016; Lei et al. Reference Lei, Lu, Liu and Wang2022), 0.55 (Hu et al. Reference Hu, Lan, Wei and Chen2019; Lan et al. Reference Lan, Hu, Yang, Wu and Chen2020) and 0.39 (Singh et al. Reference Singh, Scholl, Brinkmann, Michiel, Scheel, Herminghaus and Seemann2017). 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
$_2$) 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. Reference Ruprecht, Pini, Falta, Benson and Murdoch2014), 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
$_2$ 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. (Reference Moura, Måløy, Flekkøy and Toussaint2020) 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^\circ$, covering the wettability range of porous media from water-wet to CO
$_2$-wet. Extreme water-wet (
$\theta <45^\circ$) or CO
$_2$-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. Reference Zhao, MacMinn and Juanes2016) (see Liu et al. (Reference Liu, Berg, Ju, Wei, Kou and Cai2022) 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
$_2$ 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^\circ$. During the first CO
$_2$ injection, i.e. drainage process, the CO
$_2$ (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
$_2$ 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.
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.
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 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. Reference Ahn, Kim, Lee and Wang2020). 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).
(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
$_2$ saturation 
$S_\mathrm {CO_2}$ as a function of cycle number. The blue-dotted line with circles marks the residual CO
$_2$ saturation 
$S_{r,\mathrm {CO}_2}$ after water injection. (c) Number of CO
$_2$ ganglia. (d) Normalised average CO
$_2$ 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.
Figure 3(b) plots the CO
$_2$ saturation 
$S_\mathrm {CO_2}$ during cyclic injection. Red and blue lines, with water contact angle being 
$60^\circ$, correspond to drainage and imbibition processes, respectively. Clearly, hysteretic behaviour in 
$S_\mathrm {CO_2}$ is observed with increase/decrease in 
$S_\mathrm {CO_2}$ during drainage/imbibition. One parameter of significant importance in the context of carbon geosequestration is the residual CO
$_2$ saturation 
$S_{r,\mathrm {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,\mathrm {CO}_2}$ with number of injection cycle is observed, consistent with existing experimental and numerical observations (Edlmann et al. Reference Edlmann, Hinchliffe, Heinemann, Johnson, Ennis-King and McDermott2019; Ahn et al. Reference Ahn, Kim, Lee and Wang2020; Herring et al. Reference Herring, Sun, Armstrong, Li, McClure and Saadatfar2021). 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 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
$_2$ connecting isolated CO
$_2$ 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
$_2$ dissolution into water as a result of greater area-to-volume ratio. Note that the dissolution of CO
$_2$ 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 Reference Matter and Kelemen2009). It is found that 15 injection cycles are sufficient for the displacement process to reach the equilibrium state, i.e. no further changes in the saturation hysteresis (
${\rm \Delta} S_\mathrm {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
$_2$
First, the CO
$_2$ saturation after the first CO
$_2$ injection 
$S_{r,\mathrm {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,\mathrm {CO}_2}^{N=0}$ also represents the sweep efficiency of a single-time displacement process, which has been investigated extensively (Trojer et al. Reference Trojer, Szulczewski and Juanes2015; Jung et al. Reference Jung, Brinkmann, Seemann, Hiller, Sanchez de La Lama and Herminghaus2016; Zhao et al. Reference Zhao, MacMinn and Juanes2016; Singh et al. Reference Singh, Scholl, Brinkmann, Michiel, Scheel, Herminghaus and Seemann2017; Wang et al. Reference Wang, Chauhan, Pereira and Gan2019). The increase in 
$S_{r,\mathrm {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 Reference Cieplak and Robbins1988; Holtzman & Segre Reference Holtzman and Segre2015). Figure 4(a) shows that as the water contact angle increases, i.e. invading phase being more wetting to the porous medium, 
$S_{r,\mathrm {CO}_2}^{N=0}$ increases, consistent with past observations (Jung et al. Reference Jung, Brinkmann, Seemann, Hiller, Sanchez de La Lama and Herminghaus2016; Wang et al. Reference Wang, Chauhan, Pereira and Gan2019, Reference Wang, Pereira, Sauret, Aryana, Shi and Gan2022a). Particularly, as the contact angle increases to 
$\theta =135^\circ$, i.e. in CO
$_2$-wet porous media, a compact CO
$_2$ invasion with almost 100 % CO
$_2$ saturation is observed.
(a) CO
$_2$ saturation as a function of contact angle after: (i) the first CO
$_2$ injection, 
$S_{r,\mathrm {CO}_2}^{N=0}$ (red); (ii) the first water injection, 
$S_{r,\mathrm {CO}_2}^{N=1}$ (blue); and (iii) 15th water injection, 
$S_{r,\mathrm {CO}_2}^{N=15}$ (black). Error bars represent the standard deviation. (b) Difference in saturation between 
$S_{r,\mathrm {CO}_2}^{N=0}$ (red) and 
$S_{r,\mathrm {CO}_2}^{N=15}$, which quantifies the ultimate CO
$_2$ saturation reduction due to cyclic injection.
Then, we plot the residual CO
$_2$ saturation after the first water injection process, 
$S_{r,\mathrm {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
$_2$ saturation can be observed after the water flooding for all wetting conditions, with more significant decrease in CO
$_2$ for greater 
$\theta$. After 15 cycles, the residual CO
$_2$ saturation 
$S_{r,\mathrm {CO}_2}^{N=15}$, which represents the maximum CO
$_2$ 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
$_2$ saturation can be seen, revealing the favourable effect of cyclic injection for enhancing capillary trapping of CO
$_2$. Further, the results also indicate that at an equilibrium state, either water-wet or CO
$_2$-wet porous media are better at trapping CO
$_2$ 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 Reference AlRatrout, Blunt and Bijeljic2018). It is found that the boundary effect enhances CO
$_2$ trapping near the solid walls for CO
$_2$-wet porous media and suppresses CO
$_2$ 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
$_2$ saturation under different wetting conditions is quantified by the difference between 
$S_{r,\mathrm {CO}_2}^{N=0}$ and 
$S_{r,\mathrm {CO}_2}^{N=15}$ in figure 4(b), where an increasing trend can be observed. The insets show the CO
$_2$ (white) and water (blue) distribution after the first CO
$_2$ injection (left) and after 15 injection cycles (right) for a water-wet case (
$\theta =45^\circ$) and a CO
$_2$-wet case (
$\theta =135^\circ$), respectively. Under the water-wet condition, the initial injection of CO
$_2$ corresponds to a drainage process, where the invasion morphology falls into the capillary fingering regime described by the theory of invasion percolation (Wilkinson & Willemsen Reference Wilkinson and Willemsen1983; Lenormand & Zarcone Reference Lenormand and Zarcone1985). The resulting CO
$_2$ 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
$_2$ 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 Reference Cieplak and Robbins1988; Holtzman & Segre Reference Holtzman and Segre2015), resulting in a displacement efficiency close to 1. Nevertheless, a significant amount of CO
$_2$ can be mobilised during subsequent water flooding processes, dramatically reducing the residual CO
$_2$ saturation in the porous media, which is associated with the large value of (
$S_{r,\mathrm {CO}_2}^{N=0}-S_{r,\mathrm {CO}_2}^{N=15}$). In summary, although a CO
$_2$-wet porous medium can have larger capacity in storing CO
$_2$ during the first CO
$_2$ injection, a significant amount of CO
$_2$ 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
$_2$ saturation hysteresis during the drainage–imbibition cycles is mainly due to the increase in the residual CO
$_2$ saturation after water flooding, implying that a significant amount of CO
$_2$ is broken up by injected water into smaller blobs during the imbibition process, which leads to more stable capillary trapping. For CO
$_2$-wet porous media (
$\theta =135^\circ$), in contrast, the decrease in saturation hysteresis results from the decrease in CO
$_2$ saturation after CO
$_2$ injection. This is associated with an increasing amount of water being stably trapped within the porous media, reducing the pore space available to CO
$_2$. For neutral-wet porous media, the saturation hysteresis decreases during both water injection and CO
$_2$ 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 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
$_2$ in this case have the same contact angle 
$\theta =90^\circ$.
Capillary pressure versus CO
$_2$ 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
$_2$ and water injections. Dashed arrows as well as the increasing colour intensity indicate the increasing cycle number.
To evaluate how the saturation hysteresis 
${\rm \Delta} S_\mathrm {CO_2}$ evolves during cyclic injection, figure 6(a) plots 
${\rm \Delta} S_\mathrm {CO_2}$ against the number of injection cycles with the experimental data from Ahn et al. (Reference Ahn, Kim, Lee and Wang2020) added for comparison. It is found that the data can be described by the following exponential decay:
\begin{equation} {\rm \Delta} S_\mathrm{CO_2}=(S^{N=0}_{r,\mathrm{CO}_2}- {\rm \Delta} S^*_\mathrm{CO_2})\cdot a^{{-}N}+{\rm \Delta} S^*_\mathrm{CO_2},\end{equation}
with 
$N$ the number of cycles, 
$a$ a fitting parameter, 
$S^{N=0}_{r,\mathrm {CO}_2}$ the initial CO
$_2$ saturation after the first CO
$_2$ injection and 
${\rm \Delta} S^*_\mathrm {CO_2}={\rm \Delta} S^{N=15}_\mathrm {CO_2}$ 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 
${\rm \Delta} S_\mathrm {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.
We then plot the initial and residual CO
$_2$ saturation over the course of 15 cyclic injections in figure 7 to compare with the classic trapping model by Land (Reference Land1968). Data from two recent experimental works on cyclic injection are also added (Ahn et al. Reference Ahn, Kim, Lee and Wang2020; Herring et al. Reference Herring, Sun, Armstrong, Li, McClure and Saadatfar2021). The direction of increasing cycles are marked by dashed arrows. The movement in the initial-residual CO
$_2$ 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 Reference Land1968) (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
$_2$ 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.
Evolution of initial and residual CO
$_2$ saturation during cyclic injections for different contact angles. Dashed arrows indicate the direction of increasing cycles for 
$\theta =\{45^\circ,90^\circ,135^\circ \}$. Black-solid curves correspond to the Land model with parameter 
$C=\{0.4,2\}$.
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 
${\rm \Delta} S^*_\mathrm {CO_2}={\rm \Delta} S^{N=15}_\mathrm {CO_2}$, 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. Reference Holtzman, Dentz, Planet and Ortín2020). The terminal saturation hysteresis 
${\rm \Delta} S^*_\mathrm {CO_2}$ as a function of contact angle is plotted in figure 8. The 
${\rm \Delta} S^*_\mathrm {CO_2}$ also represents the converged active area, which is shown in the insets of figure 8 for sample cases with 
$\theta =\{45^\circ,90^\circ \,135^\circ \}$. It can be seen that 
${\rm \Delta} S^*_\mathrm {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 
${\rm \Delta} S^*_\mathrm {CO_2}$ can be explained by the symmetrical 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 Reference Hashemi, Blunt and Hajibeygi2021).
Terminal saturation hysteresis 
${\rm \Delta} S_\mathrm {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 \}$.
To interpret the hysteresis during immiscible fluid displacement in porous media, a discrete-domain model was recently proposed by Cueto-Felgueroso & Juanes (Reference Cueto-Felgueroso and Juanes2016) 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 Reference Blunt2017). 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. Reference Zarikos, Terzis, Hassanizadeh and Weigand2018) and mechanisms such as droplet fragmentation can occur under high viscous stress (Pak et al. Reference Pak, Butler, Geiger, van Dijke and Sorbie2015). 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 Reference Chiquet, Broseta and Thibeau2007; Kim et al. Reference Kim, Wan, Kneafsey and Tokunaga2012; Wan, Kim & Tokunaga Reference Wan, Kim and Tokunaga2014), 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. Reference Zhao, Kang, Yao, Viswanathan, Pawar, Zhang and Sun2018; Jahanbakhsh, Shahrokhi & Maroto-Valer Reference Jahanbakhsh, Shahrokhi and Maroto-Valer2021; Irannezhad et al. Reference Irannezhad, Primkulov, Juanes and Zhao2023), grain shape (Rokhforouz & Akhlaghi Amiri Reference Rokhforouz and Akhlaghi Amiri2019; Wang et al. Reference Wang, Pereira and Gan2021) and surface roughness (Tanino, Ibekwe & Pokrajac Reference Tanino, Ibekwe and Pokrajac2020; Wang, Pereira & Gan Reference Wang, Pereira and Gan2020; Zulfiqar et al. Reference Zulfiqar, Vogel, Ding, Golmohammadi, Küchler, Reuter and Geistlinger2020), 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
$_2$ saturation after water injection during cyclic injection, i.e. the increased stability of trapped CO
$_2$, 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
$_2$ saturation during cyclic injection. Despite that higher CO
$_2$ saturation can be achieved in CO
$_2$-wet porous media after the first CO
$_2$ injection, the majority of CO
$_2$ 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
$_2$ 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
$_2$-wet porous media, less CO
$_2$ is trapped within the neutral-wet ones at equilibrium state. In addition, the terminal CO
$_2$ saturation hysteresis, quantified by the saturation change between water/CO
$_2$ injections after sufficient cycles, is greater in neutral-wet porous media than either water-wet or CO
$_2$-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.
Appendix A. Boundary effect
To evaluate the sidewall (longer side) boundary effect on the CO
$_2$ distribution, the CO
$_2$ distribution after 15 injection cycles along the horizontal direction is calculated by summing all the CO
$_2$ 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 
$\theta =45^\circ$. 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
$_2$ occupation of pore space is less favourable along the boundary, whereas for CO
$_2$ wet porous media, more CO
$_2$ is distributed along the boundary walls. To exclude the boundary wall effect, the residual CO
$_2$ 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
$_2$ can be observed, implying the boundary effect should not influence the conclusions derived from the main text.
Boundary effects on the CO
$_2$ distribution. (a) CO
$_2$ distribution along the normalised domain width for contact angles 
$\theta =\{45^\circ \,90^\circ,135^\circ \}$. The insets show the raw data (periodically fluctuating curves) and the moving average with span length of one pore size. (b) CO
$_2$ saturation at different injection cycles (same notation as figure 4a) calculated from the area between 0.2 and 0.8 normalised domain width.
Appendix B. Pressure signal sampling and conversion
An example of the raw capillary pressure signal with 
$\theta =90^\circ$ during the first CO
$_2$ 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^*_c$ 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.
(a) An example of the raw pressure signal for 
$\theta =90^\circ$ during CO
$_2$ injection. (b) Magnified region of the first 300 computation steps. Red circles mark the sampled pressure with a positive increment.
Open access
