The impact of drainage displacement patterns and Haines jumps on CO2 storage efficiency

Injection of CO2 deep underground into porous rocks, such as saline aquifers, appears to be a promising tool for reducing CO2 emissions and the consequent climate change. During this process CO2 displaces brine from individual pores and the sequence in which this happens determines the efficiency with which the rock is filled with CO2 at the large scale. At the pore scale, displacements are controlled by the balance of capillary, viscous and inertial forces. We simulate this process by a numerical technique, multi-GPU Lattice Boltzmann, using X-ray images of the rock pores. The simulations show the three types of fluid displacement patterns, at the larger scale, that have been previously observed in both experiments and simulations: viscous fingering, capillary fingering and stable displacement. Here we examine the impact of the patterns on storage efficiency and then focus on slow flows, where displacements at the pore scale typically happen by sudden jumps in the position of the interface between brine and CO2, Haines jumps. During these jumps, the fluid in surrounding pores can rearrange in a way that prevent later displacements in nearby pores, potentially reducing the efficiency with which the CO2 fills the total available volume in the rock.

Multiphase flow is ubiquitous in nature, as well as a plethora of industrial processes. Examples include geological sequestration of CO 2 , enhanced oil recovery (EOR), water infiltration into soil, soil remediation through the removal of liquid pollutants etc. Undoubtedly, there is an inherent complexity in investigating these processes, introduced by the variety of the flow regimes, due to the interplay of the fluids' affinity to the solid surfaces (wettability) 1,2 and the complex geometry, where the fluid flow takes place, e.g. permeable media. Compact or non-compact fluid front displacement affects the displacement efficiency, defined as the fraction of the defending fluid that is displaced from the porous rock. For example the non-compact fluid front due to fingering instabilities decreases the displacement efficiency for both EOR and CO 2 sequestration. Quantifying the displacement patterns is therefore essential for optimizing subsurface processes. A combination of experiments and numerical simulations can elucidate the role of the aspects affecting multiphase flow and help construct upscaled models needed to understand the processes at larger length scales.
Our focus here is on the immiscible displacement of a wetting fluid in permeable media, termed as drainage. Extensive work in the literature addressed the above problem. According to the pioneering work of Lenormand et al. 3 in micromodels, the competition of viscous and capillary forces leads to the basic drainage displacement patterns, namely i) viscous fingering, ii) capillary fingering and iii) stable displacement. These patterns can be mapped on a phase diagram with axes the capillary number Ca and the viscosity ratio M of the fluids. Since then experimental work in fabricated micromodels [4][5][6][7][8] and numerical investigations 9-12 examined the impact of drainage displacement patterns and their domains of validity. Systematic experimental studies covered the impact of wetting conditions 13,14 , pore-scale heterogeneity 15 , as well as the phase of the injected CO 2 16 . Attention was also given on the crossover between the domains [4][5][6][7]17,18 . In micromodels 7 and in rough fractures 18 a decrease was observed in the displacement efficiency at the transition zone from viscous fingering to capillary fingering. This trend however hasn't been observed in three dimensional porous rocks 10 .
Advances in synchrotron-based X-ray computed microtomography enabled the imaging of pore-scale displacement events in porous rocks in real time without disturbing the flow 19,20 . This means that pressure gradients and the viscocapillary balance is maintained during imaging. During slow flow the pore-scale displacement process proceeds not in a smooth continuous way, but as a sequence of sharp interfacial jumps or burst events, called Haines jumps [19][20][21][22] . This is reflected by the observed sharp pressure drop in the pressure data (referred as rheons 23,24 ). Even though the average flow is at low Reynolds number, inertial effects become important over a transient amount of time during the jump events, with their experimentally observed time scale being around 1-10 ms 25,26 . Haines jumps are associated with both drainage and imbibition dynamics 21,27 , as i) at the draining site the non-wetting fluid passes through a narrow throat displacing the wetting fluid in the wider pore body, while ii) imbibition takes place in surrounding throats with the wetting fluid displacing the injected non-wetting phase. This leads to fluid rearrangement during the jumps, which supplies a significant fraction of the the necessary fluid volume for draining a pore body 19 .
Here we investigate fluid-fluid displacement patterns in a range of fluid flow regimes, by varying the capillary number and the viscosity ratios of the fluids, with the aim of better understanding the implications of these flow regimes on CO 2 geological sequestration. Supercritical CO 2 (non-wetting) is being injected deep in geological formations to displace the resident fluid (e.g. brine or oil -wetting fluid) in the pore matrix and to be safely stored there over long times. Understanding the displacement processes taking place at the pore-scale is essential in maximizing the displacement efficiency and it is of paramount importance as CO 2 geological sequestration appears to be an important tool for combating global warming.
We give particular attention to the low Ca flow regime, characterised by Haines jumps, and the impact of the associated fluid redistribution on the displacement process. This aspect has not been thoroughly investigated so far. Experimentally it is not easy to identify the imbibition sites and quantify the degree of fluid redistribution 19,20 ; rather this can be assessed from the comparison of the pore filling rates and the feed flow rate of the pump 19 . Numerically, a full Navier Stokes solver is needed to include inertial effects during the jumps. Tsuji et al. 10 investigated drainage in Berea sandstone using the color-gradient lattice Boltzmann method 28 and applying a body force (pressure gradient) to drive the fluid flow. However, this can be a limiting factor in assessing the displacement efficiency in the low Ca flow regime, as decreasing the body force to decrease Ca may lead to a pressure difference not sufficient to overcome the capillary entry pressure. In this case no flow is observed and the non-wetting phase is only capable to reach up to a given injection depth in the porous rock 9,10,29 . To the best of our knowledge only Yamabe et al. 12 report numerical work on the impact of Haines jumps on displacement efficiency. They investigated drainage in a synthetic granular rock model, in relation to CO 2 geological sequestration, and made a distinction between backward and forward Haines jumps. The authors showed that forward flowing events cause a significant drop in CO 2 saturation, whereas backward flowing events cause an increase of the CO 2 saturation per injection depth. However, the impact of these pore scale events on CO 2 displacement efficiency is not clear and believe that this investigation should be extended in several ways: 1. The domain size used by Yamabe et al. 12 is 100 3 voxels, corresponding to a physical sample size of 1 mm 3 . It is not clear whether this domain size and sample resolution used are sufficient. To investigate this point, we extend the domain size to 700 3 voxels using Graphics Processing Units (GPUs) to accelerate the computations and handle the large numerical load. We also increase the sample resolution, so that the domain size corresponds to a physical sample size of 32 mm 3 . 2. Yamabe et al. 12 drive the fluid flow using a body force. Although the authors were not explicit as to how they treat the inlet/outlet boundaries, we suspect that they apply Periodic Boundary Conditions in the direction of the flow. Here we extend this investigation by using a constant injection flow rate to drive the fluid flow (velocity boundary conditions). 3. Given the small domain used by Yamabe et al. 12 , in combination with the use of body forcing and Periodic Boundary Conditions, we believe that their results could be affected by boundary effects. Haines jumps at one end may affect the fluid flow at the inlet. Here we investigate this point by using a larger domain size in combination with constant velocity boundary conditions. 4. Yamabe et al. 12 do not report investigations of possible fluid redistribution during Haines jumps. This, however, is a very important aspect of the Haines jump phenomenon 19 . In this paper, we present computational results of extensive fluid redistribution associated with Haines jumps.
Here we demonstrate that numerically we capture the main features of the Haines jumps: i) sharp increase in the non-wetting phase velocity, ii) sharp pressure drop and iii) extensive fluid redistribution. Simulations reveal that Haines jumps can potentially decrease the displacement efficiency of the injected phase, irrespective of the type of the jump, labeled as backward or forward jump event. Irreversible fluid redistribution during the events is essentially responsible for the decrease in the displacement efficiency, as wetting phase, flowing in the direction opposite to the mean direction of the flow, can block the access to other regions of the pore space due to subsequent events affecting the displacement pathways.

Results and Discussion
We directly solve the hydrodynamic equations of motion on a three dimensional geometry reconstructed from micro-CT images of Ketton limestone 30 , see Fig. 1, using multi-GPU free energy lattice Boltzmann (LB) simulations 31-33 . The simulation system size is 700 3 lattice units (l.u) at a resolution of 4.52 μm per l.u. (physical system size 31.6 mm 3 ). Simulations cover a range of viscosity ratios M = η nw /η w and capillary numbers Ca = η nw u nw /γ, where γ, η i , u i (i = w, nw) are the surface tension, viscosity of the fluids and average fluid velocity respectively. The subscripts w, nw denote the wetting and non-wetting phases. We apply a constant injection flow rate at the inlet/outlet boundaries, i.e. the fluid flow is driven by applying velocity boundary conditions by adopting the approach proposed in 34 to two-phase flow. The numerical scheme proved to be stable over long simulation times and has been already validated in our work on Haines jumps in simplified micromodels 27 . For each viscosity ratio case we fix η i (i = w, nw) and γ, while in order to decrease M we increase the wetting phase viscosity η w . The capillary number varies in the ). Typical values for the ratio of viscous to capillary forces at the pore scale are in the range of 10 −10 -10 −3 , depending on the distance from the injection point in the well bore 35 . We note here that we keep the Ohnesorge number fixed ) in all simulations, except for two simulations at the highest Ca for logM = −1 and −2, where for numerical stability reasons we had a slightly higher Oh. The Ohnesorge number, which describes the relative importance of viscous forces to inertial and surface tension forces, is fixed for a given experimental system and can be also expressed as Oh 2 = Ca/Re. For a brine -CO 2 system, Oh is in the range of 10 −3 to 10 −2 , depending on the characteristic length scale of the system L s , as well as the fluids' properties 27,36 . L s is taken to be the average pore throat radius, as this controls the pressure at which pores drain. Finally, the equilibrium contact angle is set to θ eq = 40°, consistent with contact angle measurements in Ketton at reservoir conditions for a supercritical CO 2brine system 37 .

Drainage displacement patterns.
A qualitative demonstration of the displacement patterns in the (Ca, M) phase diagram is shown in Fig. 2, where only the non-wetting phase is shown when it reached the outlet reservoir. We verify the existence of the three typical fluid displacement patterns, namely viscous fingering (high Ca, < M log 0), capillary fingering (low Ca) and stable displacement (high Ca, > M log 0) 3 , with the main focus given on the first two. Our aim is to examine the impact of the displacement patterns on the non-wetting phase saturation: i) when the non-wetting phase reaches the outlet (breakthrough), S br nw , as well as ii) the maximum non-wetting phase saturation, S max nw , obtained if the injection continues beyond the breakthrough point. The latter determines the maximum achievable displacement efficiency of the injected phase, which in real porous media can reflect encountering a fracture or a higher permeability zone. Fig. 3(a1-c1) shows the average value for the three components of the velocity for the wetting/non-wetting phase for simulations with = M log 0 and varying Ca; the right panel, Fig. 3(a2-c2), shows the corresponding average magnitude of the velocity and the inlet/outlet pressure difference, ΔP = P inlet − P outlet . The mean direction of the flow is in the x-direction. As Ca decreases a distinct change in the flow field marks the transition to the capillary fingering regime. Sharp interfacial jumps, indicative of Haines jumps, become profound and lead to a significant increase of the fluids velocity, see Fig. 3(c). The occurrence of Haines jumps is also evident from the abrupt pressure drop observed in Fig. 3(c2), which coincides with the interfacial jumps.
Non-wetting phase saturation Vs Drainage displacement patterns. Examining the non-wetting phase saturation, S nw , as a function of the frontal position, see Fig. 4(a,b), provides useful information about the type of dispacement. The frontal position is defined as the distance from the inlet of the most advanced tip of the non-wetting phase. High Ca results in a linear increase of S nw with the frontal position of the injected non-wetting phase. As Ca decreases (transition towards capillary fingering) a distinctive step-like structure emerges. Similar behavior is observed for < M log 0 ( Fig. 4(b)) which favors viscous fingering at high Ca, with the transition from viscous to capillary fingering occuring at a lower Ca. The step-like structure observed for capillary fingering displacement is due to the sequential forward and backward Haines jump events. A forward Haines jump (or sequence of events) will increase the frontal position significantly, while S nw under the time scale of the jumps will not increase significantly. Rather what is happening is that the neighbouring regions will provide the non-wetting phase needed for filling a wider pore space through fluid rearrangement. When the interface reaches a point where the capillary entry pressure for moving forward is higher than the other available displacement pathways, for example the location indicated with the red arrow in Fig. 4(a), then jump events will be observed in the y and This leads to significant increase of S nw (Fig. 4(a) 1 → 2). The capillary fingering displacement pattern leads to loops of the non-wetting phase towards the inlet and enhances the connectivity of the non-wetting phase, as evident in Fig. 4(c) from the reduction in the Euler characteristic χ nw 38,39 . Visually this is also shown in the subfigures with different colours denoting the disconnected clusters of the non-wetting phase. We measure χ nw using a robust tool in Matlab 40 and normalize it with the total volume to describe it as a density in physical units [mm −3 ].
If injection continues, beyond the time the non-wetting phase reached the outlet (t = t br ), then S nw may increase further (up to maximum achievable S max nw ) depending on the inlet/outlet pressure difference ΔP. At high Ca the injected phase will continue to displace the defending fluid due to the high ΔP (case Ca = 3.1 × 10 −4 in Fig. 5(a)), while for low Ca no further increase in S nw will be observed (case Ca = 3.8 × 10 −5 in Fig. 5(a)), as the fluid will flow through the formed displacement pathways. Hence, once capillary fingering dictates the displacement process no significant change in S nw is achieved beyond S br nw (  S S br max nw nw ). This is evident from results shown in Fig. 5. Moreover, capillary fingering displacement commences at a higher Ca for = M log 0 (Ca ~ 4 × 10 −5 ) than for < M log 0 (Ca ~ 9 × 10 −6 for = − M log 1), as expected. Here we assume capillary fingering displacement to be commencing as . Irrespective of the viscosity ratio though, the limiting value for S max nw seems to remain the same ( for results in the capillary fingering regime. In order to minimize the impact of capillary end effects on S max nw , we also examine S nw in the first subvolume of the rock, indicated with the red arrow in Fig. 4(a) (x ≤ 1.9 mm). This is shown in the inset of Fig. 5(b) and demonstrates the same qualitative behavior.
Our results agree with the ones reported in the literature 3,6,7,12 , i.e. the maximum achievable displacement efficiency decreases with decreasing Ca and as we move from stable displacement to capillary fingering and viscous fingering. Another important observation is that the saturation at breakthrough S br nw (filled symbols) generally increases with decreasing Ca (see Fig. 5(b)), in agreement with numerical results by Tsuji et al. 10 , who investigated drainage in Berea sandstone using a color gradient LB approach 41 . Interestingly though, this is not the case for = M log 0 as Ca decreases from Ca = 3.1 × 10 −4 (Fig. 2I) to 3.8 × 10 −5 (Fig. 2J). On the contrary a significant decrease in S br nw is observed (~−17% in the whole domain and ~−12% in the first subvolume region, x ≤ 1.9 mm). Intuitively it would be reasonable to expect that S br nw should increase due to capillary fingering controlling the displacement process. What comes into play, going from Ca = 3.1 × 10 −4 (Fig. 2I) to 3.8 × 10 −5 (Fig. 2J), is the onset of Haines jumps, which becomes profound from the flow field and pressure drop shown in Fig. 3(b,c). On the other hand a decrease is observed in the inlet/outlet pressure difference ΔP, as shown in Fig. 6(a), that could justify the decrease in S br nw . Hence, the question to be addressed is whether this reduction is due to: i) the onset of Haines jumps or ii) the lower pressure difference driving the fluid flow for the lowest Ca case (Fig. 6(a)), which makes it difficult for the non-wetting phase to access regions of the pore matrix due to their higher capillary entry pressure. This will be addressed in the next section.  (Fig. 2J), Ca = 1.0 × 10 −5 (Fig. 2K) and ii) logM = −1: Ca = 8.5 × 10 −6 (Fig. 2G). Particular focus is given on the fluid redistribution associated with the events and how this can affect the displacement process. As demonstrated experimentally 19,21 and numerically 27,42 , Haines jumps are associated with both drainage and imbibition dynamics, as wetting phase imbibes surrounding regions of the jump, displacing the non-wetting phase and providing so a significant volume of non-wetting phase to the draining pores. To assess this quantitatively, we identify the pore filling events from the pressure signal, see Fig. 6(b). Using the fluids' configuration at the beginning (t start ) and the end (t end ) of each event we quantify the degree of fluid redistribution R % , as we can measure directly the draining volume V d , as well as the non-wetting phase volume originating from the imbibition sites V imb . The degree of fluid redistribution, defined as R % = V imb /V d , is shown in Fig. 6(c). As the injection flow rate (Ca) decreases for logM = 0 (J to K), R % increases from R % ~ 39 ± 12% up to even R % ~ 80% (63 ± 12%). Decreasing the viscosity ratio to logM = −1 (K to G by increasing η w and keeping u inj nw fixed) results in a decrease in R % (50% ± 12%), reflecting the increased viscous dissipation rate in the wetting phase in G. Nevertheless results verify that extensive fluid redistribution takes place, in line with the experimental findings that the required draining volume can not be supported by the externally imposed flow rate in the time scale of the jump 19,21 . Berg et al. 19 report R % ~ 99% in Berea sandstone at even lower Ca ~ 10 −8 . Moreover, as expected, a  (Fig. 2G) can be visualized in Fig. 6(d). The region occupied by the non-wetting phase that remains unchanged during the event is  A rough estimate for the time scale of the interfacial jumps, Δt j = t end − t start , considering the capillary forces due to differences in the interfacial curvature and the mass of the accelerated fluid 25 3 l.u, where β > 1 is a dimensionless constant that depends on the pore shape. This analysis neglects the viscous resistance to the flow; hence it is justifiable to obtain a higher time scale in the simulations. Results are an order of magnitude higher for logM = 0 (Δt j = 4.5 ± 2.4 × 10 4 − Ca = 3.8 × 10 −5 , Δt j = 6.5 ± 4.9 × 10 4 − Ca = 1.0 × 10 −5 ), while as the viscous resistance in the wetting phase increases further for logM = −1, Ca = 8.5 × 10 −6 , we obtain Δt j = 10.5 ± 5.2 × 10 4 .
From an energy point of view, during the drainage process surface energy is stored in the system due to the externally performed work of pressure, However, due to viscous dissipation, the change in surface energy ΔF surf is less than the work done on the system 44 . This is shown in Fig. 7(a) where we plot the efficiency of the conversion of work of drainage (pressure-volume work) to surface energy, E d = ΔF surf /W p , as a function of Ca. At the highest Ca a significant amount of energy is also converted to kinetic energy E k ≥ 0. As the injection flow rate drops 3 orders in magnitude ( → − − u : 10 10 inj nw 3 6 ), E k becomes negligible and the energy input into the system is converted mainly to surface energy and dissipated. E d increases as Ca (injection flow rate) decreases due to lower viscous dissipation rate. Seth and Morrow 44 estimate E d as a function of wetting phase saturation S w in centriguge experiments, but do not quantify the importance of viscous forces to surface tension (Ca). They report E d ~ 55% for Berea sandstone at S w ~ 0.6, which increases to E d ~ 90% for sphere packs. Our results at this value of S w reveal E d ~ 68% for the lowest Ca simulations. In the capillary fingering regime, energy stored in the menisci and the fluid columns of the non-wetting phase in the pore throats is released during the jump events, converted into kinetic energy and interfacial energy and dissipated. Hence, the system is essentially using existing energy already stored as surface energy to support these abrupt events and the fluid redistribution, while the fluids acceleration increases the energy losses due to the higher viscous dissipation. The change in surface energy is given by dF surf = γdA int + γ ws dA ws + γ ns dA ns , where dA int , dA ws and dA ns are the increments of the areas of the fluid-fluid, solid-wetting phase fluid and solid -non wetting phase fluid interfaces respectively and γ,  Fig. 2G). Inset: The draining volume, normalized by the average pore volume, as a function of pressure drop Δp j . (d) Haines jump event and the associated fluid rearrangement from a simulation at logM = −1, Ca = 8.5 × 10 −6 . The region occupied by the non-wetting phase that remains unchanged before and after the jump event is shown in yellow, while the draining pore body (V d ) and interfacial recession (V imb ) is shown in light blue and red respectively. The interfacial energy created corresponds to 56% of the total available energy (externally performed work of pressure and the energy released at the imbibition sites), while the volume of fluids being redistributed corresponds to 71% of the draining volume.
Examining the impact of Haines jumps on the displacement efficiency revealed that these events can potentially decrease the displacement efficiency if the fluid redistribution at the imbibition sites can be characterized as irreversible. This irreversibility means that after the jump the non-wetting phase never manages to displace the wetting phase back at the imbibition sites, see Fig. 8(a) with results from a simulation at logM = −1, Ca = 8.5 × 10 −6 (Fig. 2G). The wetting phase eventually becomes trapped there (t = t br ). It was observed in our previous work in micromodels 27 (see Fig. 11 therein) that reversal of the flow (imbibition) can happen even after the non-wetting phase passed the narrowest restriction point in the throat (hence exceeded the capillary entry pressure) or even escaped in the wider pore body. Therefore, it is not only the volume of the throats occupied by the non-wetting phase prior the event that matters, but most importantly the fact that the displacement pathways can be affected; this may have a bigger impact on the displacement efficiency as regions of the pore matrix become inaccessible to the non-wetting phase. Comparing the fluids' distribution prior an event and just before breakthrough, enables us to identify the total locations that can be characterized by irreversible displacement. This is illustrated in Fig. 8(b) (regions in red). The volume of these regions corresponds to a loss of ~4% of the existing pore space occupied by the non-wetting phase at the beginning of the jump event.
Another mechanism and important evidence of how Haines jumps can potentially lead to a decrease in the displacement efficiency comes from analysing an event at logM = 0, Ca = 1.0 × 10 −5 (Fig. 2K), see Fig. 9. Initially two jump events take place in the two locations shown in Fig. 9(a2) at t = t 2 . The jump event, indicated in red cycle, proceeds in a cascade like manner over multiple geometrically defined pores, causing a big reduction in capillary pressure even in throats that are a significant distance away from the jump event. The reduction in the capillary pressure in this throat below its snap-off capillary pressure causes the disconnection of the non-wetting phase and produces a long-lasting fluid configuration as it remains disconnected until the end of the simulation. Distal snap-off events may have a bigger impact on the flow than local snap-off events, and have been observed experimentally for CO 2 -brine system 43 . Here for example, due to the snap-off event, the injected phase loses access to the pore space shown in light blue in Fig. 9(b), which is filled by the non-wetting phase for = M log 0, Ca = 3.1 × 10 −4 (Fig. 2I). Therefore, such events can firstly affect the displacement pathways and potentially lead to a decrease in the displacement efficiency.
The above two mechanisms demonstrate how Haines jumps potentially affect the displacement efficiency. Irrespective of the type of Haines jumps, whether backward or forward, these abrupt events decrease the storage efficiency of the non-wetting phase, contrary to what is argued by Yamabe et al. 12 . Their findings refer to the impact of the jump events on CO 2 saturation per injection depth, similar to what is shown in Fig. 4, but do not actually examine the impact of Haines jumps on storage efficiency, nor the aspect of fluid rearrangement. Furthermore, an important remark to make here is that using a constant injection flow rate (velocity boundary conditions) to drive the fluid flow enables the study of the low Ca flow regime characterised by Haines jumps. Even as the injection flow rate decreases, pressure slowly builds up in the system until the capillary entry pressure is overcome. This would not have been possible if a body force (pressure gradient) is used, as the injected phase may not be able to reach the outlet of the domain if the overall pressure difference generated does not overcome the required capillary entry pressure along the percolation pathway 9,10,29 .

Conclusions
Drainage displacement patterns have been identified in the pioneering work of Lenormand et al. 3 , as capillary fingering (low Ca), viscous fingering (high Ca, < M log 0) and stable displacement (high Ca, > M log 0). These distinctively different displacement flow regimes can affect significantly the displacement efficiency, which is defined as the fraction of the defending wetting fluid that has been displaced from the pore matrix when the injected non-wetting phase reached the outlet of the domain. A more important feature is the maximum achievable displacement efficiency, established by continuing the injection until saturation convergence is achieved.
Here, in order to understand these pore-scale displacement phenomena and their impact on CO 2 storage efficiency, we investigate the immiscible displacement (drainage) of a wetting fluid in a porous medium by a non-wetting fluid using multi-GPU free energy lattice Boltzmann simulations, under various capillary numbers Ca and viscosity ratios M. As a first step we first verify the existence of the three typical fluid displacement patterns, before turning our attention on the impact of these processes on storage efficiency. Our results agree with the ones reported in the literature 3,6,12 , i.e. maximum achievable displacement efficiency decreases with decreasing Ca and as we move from stable displacement to capillary fingering and viscous fingering.
Then we focus on flow at low Ca (capillary fingering regime), and the impact of Haines jumps on displacement efficiency, by decreasing the injection flow rate. During these sharp interfacial jumps pressure drops abruptly. Pressure drop coincides with the sharp increase of the non-wetting phase velocity. Extensive fluid rearrangement, which increases as Ca decreases, provides extra non-wetting phase needed for draining wider pore bodies. We identify two possible mechanisms that potentially affect the displacement process: i) irreversible imbibition displacement and ii) distal snap-off events. Both can render regions of the pore space inaccessible to the injected non-wetting phase. During the redistribution of the fluids, wetting phase flowing in the direction opposite to the mean direction of the flow displaces the non-wetting phase. If this imbibition displacement is irreversible, it can lead to wetting phase trapped due to subsequent jump events, but most importantly leads to an overall decrease in the displacement efficiency by subsequently blocking access to other regions of the pore space. Distal snap-off events can also change the displacement process, as the disconnection of the non-wetting phase may prohibit the injected phase from accessing regions of the porous matrix. From an energy point of view, the extensive fluid redistribution associated with the events and the increased viscous dissipation rate as the fluids accelerate, mean that less energy is eventually stored in the system as surface energy; hence this decreases the efficiency of the conversion of work of drainage to surface energy.
Hence, our findings have important implications in the context of geological sequestration of CO 2 , as Haines jumps become a limiting factor in the storage process. Suppressing these events, for example by decreasing surface tension through the use of surfactants, can shift the flow regime to higher capillary numbers and increase considerably the saturation of CO 2 that can be attained in the porous matrix. Eventually, it would be desirable to restore surface tension to its original value, e.g. adsorption of surfactants at solid surfaces or injection of brine to ⁎ t 0 862 end ), as well as the corresponding fluid redistribution. The region occupied by the non-wetting phase prior the jump but not after is shown in red (imbibition sites). Subsequent events don't change the occupancy of this pore space and the wetting phase becomes eventually trapped there (irreversible imbibition displacement), see configuration at breakthrough (t = t br ). (b) Total irreversible imbibition sites, comparing the fluids' distribution at t = 18.18 × 10 6 (prior a jump, t* = 0.856) and t = 21.15 × 10 6 (before breakthrough, t* = 0.996). The region occupied by the non-wetting phase that remains unchanged is shown in yellow, while the draining regions and interfacial recession sites (irreversible displacement) are shown in light blue and red respectively.
SCIEntIfIC REPORTS | (2018) 8:15561 | DOI:10.1038/s41598-018-33502-y dilute the concentration of surfactants, as high capillary forces are needed to trap CO 2 in the pore space over long time scales. The latter can be in line with water-alternating-gas injection schemes 45 and the findings of Herring et al. 38 who demonstrated that both the drainage (CO 2 injection) and imbibition (subsequent water injection) processes can be engineered in order to maximize residual trapping within the permeable medium.

Methods
Porous rock. The three dimensional geometry used for the simulations was reconstructed from micro-CT images of Ketton limestone (porosity 0.159 and permeability 9.4 Darcy) 30 . We note that the Ketton limestone has microporous grains with micro-pore sizes well below the resolution of the CT scan 46 and this is neglected in the flow simulations presented in the paper. However, given that we expect the permeabilities of the microporosity and macroporosity regions to be several orders of magnitude apart, and that we are considering capillary pressures below the entry pressure of the micro-porous regions, the proposed approach to simulate flow in the pore space that corresponds to the macroporosity only is reasonable.
Numerical Model -Free energy lattice Boltzmann method. Here we describe the numerical method we use, covering the thermodynamics, the dynamical equations of motion and finally the numerical implementation. We have chosen to use a standard free energy lattice Boltzmann approach 32 to solve the hydrodynamic equations of motion for two-phase flow directly on micro-CT images of porous media. The method belongs to a class of hydrodynamic models, called diffuse interface (phase field) models [47][48][49][50][51] , where the sharp-interface formulation is replaced with a continuous variation of an order parameter over a finite-sized interfacial region. These models introduce a diffusive mechanism at the interfacial region, which regularizes the viscous dissipation singularity 52 , and allows the contact line to slip on a solid substrate.
Thermodynamics of the fluid. The equilibrium properties of a binary fluid can be described by a Landau free energy functional 32 . The jump event indicated in red cycle proceeds in a cascadelike manner over multiple geometrically defined pores resulting in a big reduction in capillary pressure even in throats that are a significant distance away from the jump event. The reduction in the capillary pressure in this throat below its snap-off capillary pressure causes disconnection of the non-wetting phase (distal snapoff). (b) Same distal snap-off event from a different viewing angle. The region in light blue on the right panel is pore space that, at breakthrough, is occupied by the non-wetting phase for Ca = 3.1 × 10 −4 and = M log 0 (Fig. 2I).
SCIEntIfIC REPORTS | (2018) 8:15561 | DOI:10.1038/s41598-018-33502-y The first term in the integrand is the bulk free energy density given by where φ is the concentration or order parameter, ρ is the fluid mass density and c is a lattice velocity parameter. A is a constant with dimensions of energy per unit volume. The bulk free energy density f b has minima at φ eq = ±1, corresponding to the two bulk fluid phases, while the locus φ = 0 is chosen as the position of the interface. The term in ρ controls the compressibility of the fluid 53 . The gradient term, φ ∂ κ α φ ( ) 2 2 , accounts for the excess free energy due to the presence of interfaces by penalizing spatial variations of the order parameter φ. This gives rise to the interface tension γ κ = φ A 8 /9 and to the interface width ξ κ = φ A / 32 . The affinity of the fluids to solid surfaces is controlled by the final term in the free energy functional, eq. 1. Following Cahn 54 , the surface energy density is taken to be of the form f s = −hφ s , where φ s is the value of the order parameter at the surface. Minimisation of the free energy gives an equilibrium wetting boundary condition 32 The value of the parameter h (the surface excess chemical potential) is related to the equilibrium contact angle θ eq via 32 eq where α θ = arccos(sin ) eq 2 and the function sign returns the sign of its argument. This choice of the free energy leads to the chemical potential