The Impact of Wettability and Surface Roughness on Fluid Displacement and Capillary Trapping in 2‐D and 3‐D Porous Media: 1. Wettability‐Controlled Phase Transition of Trapping Efficiency in Glass Beads Packs

Fluid invasion, displacement of one fluid by another in porous media, is important in a large number of industrial and natural processes. Of special interest is the trapping of gas and oil clusters. We study the impact of wettability on fluid pattern formation and capillary trapping in three‐dimensional glass beads packs (dmean = 1 mm) during fluid invasion at capillary numbers of 10−7 using μ‐CT. The invading fluid was water, and the defending fluid was air. The contact angle of the glass beads was altered between 5° and 115° using Piranha cleaning and silanization. We analyzed the front morphology of the invading fluid, the residual gas saturation, the fluid occupation frequency of pores, and the morphology and statistics of the trapped gas clusters. We found a sharp transition (crossover) at a critical contact angle θc = 96°. Below θc the morphology of the displacement front was flat and compact caused by the strong smoothing effect of cooperative filling. Above θc the morphology of the displacement front was fractal and ramified caused by single bursts (Haines jumps). Across this dynamical phase transition the trapping efficiency changes from no trapping to maximal trapping. For θ > θc the experimental results show that invasion percolation governs the fluid displacement. Strong indicators are the universal scaling behavior of the size distribution of large clusters (relative data error εdata < 1%) and their linear surface‐volume relationship (R2 = 0.99).


Introduction
Capillary trapping of gas clusters, residual NAPL clusters, and oil ganglia plays a central role in many natural and technological processes. Gas trapping within the capillary fringe of the unsaturated soil zone during infiltration is an important oxygen-enrichment process of surface groundwater. Capillary trapping of CO 2 in underground gas reservoirs (e.g., saline aquifers and carbon capture and storage) and direct gas injection into aquifers with subsequent gas storage, for example, H 2 , CH 4 for energy storage, and O 2 for groundwater remediation, are examples of technological processes. Most of the displacement processes are slow (quasistatic), and viscous forces can be neglected, that is, capillary numbers Ca are smaller than 10 −6 .
In their pioneering work, Chatzis et al. (1982) conducted a comprehensive study on residual oil trapping using different porous media (glass beads packs [GBS], sandstones, and 2-D micromodels [MMs]). They discussed different trapping mechanisms in different porous media and developed a clear physical picture for the two main trapping mechanisms by instructive examples: (i) Bypass trapping is visualized by pore doublets, and (ii) snap-off trapping (capillary instability) is visualized by an undulating capillary with high pore-throat aspect ratio. Their experimental results demonstrated the effect of pore structure, topology, and heterogeneity on capillary trapping.
Surface roughness of the inner pore-solid interface also influences the capillary trapping efficiency. Glass beads and natural sands show a significant difference of trapped gas phase of about 15% , although both media exhibit nearly the same grain size and pore size distribution (PSD), and connectivity as measured by μ-CT. Therefore, the critical question is: What is the physical reason for the different capillary trapping efficiency of glass beads and natural sands? Precursor thick-film flow, driven by the high capillary pressure of surface cavities of rough surfaces (sand and sandstones), leads to an efficient snap-off trapping mechanism. On the other hand, bypass trapping governs the capillary trapping process in glass beads with a smooth surface.
Recent research focused on the impact of wettability on fluid pattern formation and on capillary trapping (Wang et al., 2019). In their classical work Cieplak and Robbins (1988) discussed an interesting dynamical phase transition associated with changes in local growth modes at the displacement front. This transition also displayed fluid-fluid pattern formation in dependence of the contact angle θ of the invading phase during quasistatic displacement. In the following, we denote the observed dynamical phase transition as Cieplak-Robbins transition or percolation transition. Below a critical contact angle (θ c~6 0°; 2-D estimate; Cieplak & Robbins, 1988), that is, for smooth surfaces, when the invading fluid is the wetting fluid, the fluid floods the porous medium in a compact pattern without any capillary trapping of the defending phase (nonwetting phase; Figure 1a). There appears to be a strong smoothing force analogous to surface tension at long length scales. This smoothing of the fluid interface results from cooperative filling by neighboring porescalled "overlap" displacement process. Furthermore, the smoothing is absent in the invasion percolation (IP) model, which predicts a fractal or heterogeneous interface (displacement front). Above θ c the invading fluid is nonwetting and the drainage displacement process is driven by single throat "bursts" or Haines jumps. This produces a fractal interface, which, in turn, causes capillary trapping of the defending wetting fluid (Figure 1b). The fluid-fluid pattern generation, the fractal front behavior, and the statistics of the trapped clusters obey IP. Cieplak and Robbins (1988) numerically studied this transition for horizontal flow in a 2-D model of random-sized discs at regular lattice (triangular and quadratic) positions. The statistics of the pore-scale displacement processes show that overlaps and bursts dominate the displacement process below and above θ c , respectively. The authors argued that the greater connectivity of interfaces in 3-D will cause (i) a larger critical contact angle and (ii) a stronger smoothing effect.
In previous work, we observed this transition in 3-D using a glass beads monolayer (ML) for vertical flow with gas as the defending fluid and water as the invading fluid . The ML represents a quasi-3-D porous media, where randomness results from the heterogeneous wettability of the silica surface. Figure 2 shows the case of θ~20 ± 5°< θ c . At constant flow rate (Ca = 10 −6 ) water searches for a flow path with lowest flow resistance, that is, for a gas-water interface with small curvature. The gas-water interface (denoted as Front in Figure 2) is visualized by the thick black curve. At initial time t = t 1 the front is flat, and at a random position the front has to move upwards. The white interface moves continuously with Figure 1. Front movement in z direction against gravity and bypass gas trapping. The stable menisci after Haines jumps at subsequent times (1, 2, and 3) are shown by red lines. The solid and dashed black lines show the menisci during pressure increase. The dashed lines give the last stable menisci before they "burst" or "jump" into a more stable configuration (2). (a) Water-wet sediment (contact angle θ < θ c ): Overlap of two neighboring menisci leads to cooperative filling and interface smoothing. The compact front pushes the gas bubble (yellow disc) away, that is, the chance of bypass gas trapping is dramatically reduced. (b) Intermediate-wet sediment (θ > θ c ): The two neighboring front menisci do not overlap during front advance. This separate meniscus advance causes bypass trapping of a gas bubble (3).
increasing water pressure to the unstable position z 2 (low capillary pressure) and then "bursts" (jumps) to z 3 (Haines jump). A further vertical move of the front has to overcome the high capillary pressure given by the curvature (∝1/r 1 ) of the white interface. An alternative pathway with lower capillary pressure (lower curvature) exists in the lateral direction. Cooperative filling of adjacent pores of Sphere 1 leads to a smaller curvature (∝1/r 2 ) compared to vertical displacement process. Note that this searching for pathways with lower capillary pressure in this regular ML structure is the same as for cooperative filling in random sphere packs. Figure 3a shows this lateral displacement process and the strong smoothing of the fluid interface caused by cooperative filling. As the gas is pushed away by this lateral front movement, no significant gas trapping occurs, only small tiny gas bubbles (red points) are trapped, always in direction of the front advance. Figure 3b shows the transition to the intermediate-wet sediment (θ = 70 ± 20°> θ c ) characterized by a fractal interface and strong capillary trapping. Both the cluster size distribution and the correlation length can be predicted by IP theory: (i) the cluster size distribution exhibits a universal power law behavior with an exponent τ exp = 2.06 that deviates from the theoretical 3-D value by 5% only . This good agreement indicates that the capillary trapping within the 2-D ML is governed by the 3-D critical exponent and, therefore, behaves similar to a 3-D porous media. (ii) The experimental data for the cluster heights show the significant influence of buoyancy; therefore, a natural length scale-a correlation length or cutoff length-is introduced by the buoyancy force as predicted by the IP theory (Wilkinson, 1984). The 3-D correlation length yields an upper limit for almost all of the trapped clusters. The videos of the Cieplak-Robbins transition are available online (https://www.sciencedirect.com/science/article/abs/pii/S0021979715300941-mmc1.mp4 and https://www.sciencedirect.com/science/article/abs/pii/S0021979715300941-mmc2.mp4).  There are contradicting experimental and numerical results about the impact of wettability on the capillary trapping efficiency. Singh et al. (2017), Hu et al. (2017), Jung et al. (2016), and Trojer et al. (2015) found an increasing capillary trapping efficiency. In contrast, Herring et al. (2016), Rahman et al. (2016), Chaudhary et al. (2013), andIglauer et al. (2012) found decreasing capillary trapping efficiency with increasing contact angle.
This discrepancy is presumably caused by the surface roughness of the inner pore-solid interface. Therefore, it is necessary to distinguish between the impact of wettability on sediments with smooth surface (glass beads, silica MMs, Si MMs, and GBS ML) and sediments with rough surfaces (sands, sandstones, natural rocks, and MMs produced by glass ceramics). As the Cieplak-Robbins theory applies to front advance caused by menisci (full duct flow), the predictions are valid for porous media with smooth surfaces. In this case, the relevant trapping mechanism is bypass trapping. Singh et al. (2017) studied the displacement process at Ca = 1.6 × 10 −7 in dense random packs (porosity = 0.39) of glass beads (355-425 μm) by real-time μ-CT. The invading fluid was water, and the defending fluid was oil (dodecane, interface tension = 47 mN/m, viscosity = 1.383 mPa/s, and density = 761 kg/m 3 ). The contact angle of the glass beads was altered between 60°and 165°by silanization, different cleaning procedures, and different fluids. The Cieplak-Robbins transition was observed between 80°and 110°, that is, a critical contact angle of about 95°. The crossover in fluid pattern at the displacement front also caused a crossover in capillary trapping efficiency of the defending fluid from 0% to 13%. This means that no oil was trapped for contact angles from 0-70°. The authors derived a theoretical critical contact angle of about 90°at which these bursts occur, based on simple geometric arguments assuming ordered closest cubic packing of equal-sized glass spheres. However, the glass beads are not equal sized and are distributed between 355 and 425 μm with an average grain size diameter of 383 ± 25 μm. Furthermore, the experimental porosity of 0.39 is about 35% larger than the theoretical porosity of 0.25. This leads to random packing with disordered layer structure and altering interface angles (angle between two neighboring menisci) between 120°(triangular packing) and 180°(quadratic packing; Figure 1a). We note that the critical contact angle in a stochastic porous medium is an ensemble average over many realizations. Cieplak and Robbins (1990) emphasized that a proof of this pattern transition has to show that the fluid pattern above θ c is controlled by IP with the corresponding universal exponents, for example, for the cluster size distribution. This IP proof is currently missing for packed glass beads.
The objective of this paper is to study the Cieplak-Robbins transition within packed glass beads with different degrees of wettability using μ-CT.

Critical Contact Angle for the Cieplak-Robbins Transition
It is instructive to study the influence of random packing on the critical contact angle θ c , because random packing of spheres with a size distribution R ¼ R ± 0:1R and a porosity between 0.36 and 0.37 ( Figure 4) deviates dramatically from ordered closest cubic packing (porosity = 0.25).
We discuss the impact of short-range order and porosity on θ c , based on a 2-D triangular lattice with equilateral triangles. Cooperative filling between two neighboring menisci will occur, before the single meniscus reaches the top sphere at point P1 ( Figure 5), if the two menisci with high curvature overlap at the contact point P1. The two menisci will then form a common meniscus (interface) of smaller curvature, that is, lower capillary pressure and hence higher water phase pressure. Based on these two constraints, we calculated the critical contact angle for the Cieplak-Robbins transition ( Table 1).
The results show that both the short-range order and the porosity have a significant impact on the critical contact angle. For the experimental porosity of 0.37 (last row in Table 1), θ c varies between 8°and 64°, depending on the short-range order (interface angle α). As the experimental critical contact angle is roughly 90°, this demonstrates that regular layered packing is not realized in random bead packs. Rather, an interface angle of ≤120°dominates the packing, and the short-range order is determined by the arrangement of neighboring spheres with different radii (Figure 4b). Spheres with different radii likely also have an impact on the critical contact angle ( Figure 5).

IP Above the Critical Contact Angle
In usual percolation, the approach to the percolation transition is universal. It depends only on the dimension, not on the network structure or whether one considers percolation of sites or bonds. Both bond percolation, which corresponds most closely to nw invasion (θ > θ c , nw-non-wetting), and site percolation, which corresponds most closely to w invasion (θ < θ c , w-wetting), result in fractal invasion patterns with pore-scale structure (Blunt & Scher, 1995;Stauffer & Aharony, 1994;Wilkinson, 1984).
Comparing percolation in 2-D and 3-D porous media, a fundamental dimensionality problem arises as Wilkinson (1984) noted: "Two species (fluids) cannot percolate each other in 2D-micromodels … as this is the case in 3D-systems." In other words, phase continuity for both fluids is not possible in 2-D (Adler & Brenner, 1988;Hunt & Sahimi, 2017). This limits the value of the experimental data obtained from 2-D MM experiments to test percolation theory with trapping. In 3-D there are two percolation thresholds, one for the wetting and one for the nonwetting phase. There also is a range of phase saturations where both phases percolate (not possible in 2-D). During fluid invasion, trapping takes place while the percolating cluster is growing in a statistically homogeneous and isotropic fashion. Trapping proceeds until the percolation

Water Resources Research
GEISTLINGER AND ZULFIQAR threshold of the defending fluid is reached. Blunt and Scher (1995) denoted this threshold as the "terminal point." We emphasize that most of the trapping occurs during the time where both fluids span a percolating cluster over the porous media. Hence, phase continuity must hold for both fluids. We previously tested this fundamental statement by comparing the theoretical two-threshold picture of the trapping process with the experimental thresholds, in 1-mm GBS MLs and 2-D MMs . We note that the characteristics of the percolation transition hold only near the threshold of the defending phase, that is, at the terminal point.
If capillary forces are the dominant forces at the pore scale, and viscous forces can be neglected, the flow patterns can be described by the IP theory. We estimate the ratio of viscous forces to capillary forces for a hypothetical gas bubble trapped in a single d k,max pore. The gas bubble is held by capillary forces that are determined by a d k,min throat. Both characteristic diameters can be estimated using equations (2) and (3) in Geistlinger et al. (2006): The experimental porosity of 0.37 gives an averaged coordination number ⟨Z⟩ = 8.47 that leads to ξ min = 0.31 and ξ max = 0.62. Then the capillary pressure in a pore throat can be estimated: where σ gw denotes the gas-water interfacial tension. The viscous pressure drop across d k,max is given by where q w denotes the volumetric flux of the wetting fluid, μ w is the viscosity, and k is the absolute permeability assuming a relative permeability of order 1. Therefore, the ratio of viscous to capillary forces is with the experimental contact angle θ = 99°(Exp5). This assumes an empirical relationship between permeability and grain size for random close packing of spheres of radius r k , k ¼ δ · r 2 k , and δ ≈ 2.7 × 10 −3 (Blunt & Scher, 1995). Due to the large contact angle (|cosθ| = 0.17), the conventional capillary number (Ca = μ w v w /σ gw ) gives almost the correct ratio of viscous forces to capillary forces. For θ = 5°(Exp1, water-wet GBS), the conventional capillary number overestimates the viscous forces at pore scale by 1-2 orders of magnitude. As Exp5 was conducted at a capillary number Ca ≅ 10 −7 , capillary forces dominate viscous forces by 7 orders of magnitude. Consequently, the IP theory seems suitable for describing the experimental results for θ > θ c .

Cluster Statistics
Following Stauffer and Aharony (1994, Ch. 2.2), the probability that the trapped cluster contains exactly s sites (nodes, pores; s is given in voxels) is PDF(s) = N(s)/N tot , where N(s) is the number of residual clusters of size s counted and N tot is the total number of pore space voxels. The cumulative cluster size distribution is then given by At the percolation threshold p c (or the terminal point, at which the non-wetting fluid becomes  Blunt & Scher, 1995), and for large clusters, the normalized cluster number follows a universal power law decay with the Fisher exponent τ (Fisher, 1967): Accordingly, where k is an arbitrary constant. The critical exponents for 3-D porous media that are relevant for the following discussion are as follows (Stauffer & Aharony, 1994): τ = 2.18 and ν = 0.88.

Linear Relationship Between Cluster Surface Area and Cluster Volume
To understand the linear relationship between the total gas surface area and the trapped gas volume of the nw phase, we first study the relationship between the surface area A i of a single trapped gas cluster and its volume V i . Consider a cluster of size s that forms a spherical bubble. The surface area is expected to be proportional to the 1/D root of its volume: That means, for dimension D = 3, the exponent p = 2/3 < 1. Large clusters are flat and ganglia-like ( Figure 7c). Therefore, they exhibit a large surface area and a minimal volume. Thus, for large (but finite) clusters, the following linear relationship is expected: that is, the exponent p = 1, and hence, a linear surface-volume relationship for large trapped clusters. We note that Stauffer and Aharony (1994) provided proof that the surface area is always proportional to its volume for sufficiently large clusters (see equation 44b in Stauffer & Aharony, 1994). Therefore, this linear relationship is universal, that is, independent of pore-scale properties like bond or site percolation, PSD, or coordination number.

Influence of Buoyancy and Cutoff Length
To study the effect of buoyancy on the size of trapped gas clusters, we follow Wilkinson (1984) and Blunt et al. (1992). We therefore return to the physical picture discussed above of a gas bubble trapped in a single d k,max pore (buoyancy force = Δm · g) and held by capillary forces that are determined by d k,min , (capillary force = πd k,min σ gw cosθ). In this scenario, the ratio of buoyancy forces to capillary forces (called modified Bond number e B) for Exp5 is given by where B denotes the conventional Bond number defined by the first equality in equation (11). The occupation probability is reduced by the hydrostatic gradient Δρ · g · z (Δρ-density difference of fluids; ggravitational acceleration) for water invasion. Hence, the occupation probability varies linearly with the cluster height z. Consequently, the buoyancy force will affect the most extensive clusters and truncate PDF (p c ,s). Physically, this truncation leads to a maximum correlation length ξ B in the system according to with α = ν/(1 + ν) = 0.47, where ν = 0.88 is an universal critical exponent (Stauffer & Aharony, 1994). The interpretation of the length ξ B in the imbibition case is that the length provides an upper limit of the largest trapped clusters of the nw fluid in a statistical sense. That means that most of the trapped clusters will obey the limit. Calculating the correlation length according to equation 12 yields approximately 2 mm.

Wettability Alteration
We studied six different GBS: strongly water-wet GBS, weakly water-wet GBS, three GBS of intermediate wettability, and weakly air-wet GBS. Untreated GBS show a contact angle θ = 33°and represent, therefore, the weak water-wet GBS. Strong water-wet GBS (θ = 5°) were produced by treating GBS with Piranha solution (H 2 SO 4 :H 2 O 2 = 3:1), rinsed with excess deionized water, ultrasonically cleaned, and dried in an oven. The three intermediate-wet GBS (θ = 81°, θ = 88°, and θ = 99°) were produced by a self-assembled ML of dichlorodimethylsilane (DCDMS) that was coated on Piranha-cleaned GBS surfaces (Dickson et al., 2006). DCDMS was prepared in pure cyclohexane, and the glass beads/plates were submerged into this solution. Supernatant was replaced by pure cyclohexane, followed by washing with absolute ethanol and drying. DCDMS concentrations between 0.1 and 10 mol/m 3 alter the contact angles of the glass surface between 75°and 105°. Weakly air-wet GBS (θ = 115°) were achieved by coating Piranha-cleaned GBS with a self-assembled ML of octadecyltrichlorosilane (OTS), prepared in bicyclohexyl (Lessel et al., 2015). The glass beads/plates were posttreated with chloroform and dried. Silane concentrations and adsorbed water content on the silica surfaces are always crucial for achieving a smooth continuous ML (Zhang et al., 2008).
We used the Analog Glass Plate Method (Herring et al., 2016) to control the contact angle; that is, we measured the contact angle on a flat silica glass substrate (exhibiting identical chemical composition as silica GBS: SiO 2 = 72%, Na 2 O = 14%, CaO = 7.5%, and K 2 O = 1%). Therefore, we applied each surface treatment method to the glass plates and to the glass beads. We thereby assumed that the air-water contact angle is similar for glass beads and flat glass substrates. The air-water contact angle was measured using Drop Shape Analyzer 100 (KRÜSS, Germany).

Experimental Setup and Methodology
For the fluid invasion experiment, a Plexiglas column of 200-mm length and 19.1-mm inner diameter was packed with the GBS (Figure 6). Small vibrations were applied to the column during the continuous filling process to guarantee close and reproducible random packing, indicated by a highly ordered structure of spheres near the wall ( Figure S1). For all experiments, close random packing with approximately identical PSD (0.273 ± 0.005 mm) was achieved. The PSD was derived by the first Minkowski function of the pore space . Note that the standard deviation is smaller than 2%, which indicates a good reproducibility of the packing procedure. To avoid grain rearrangement during fluid invasion, lead spheres (diameter = 3.1 mm and density = 11.35 × 10 3 kg/m 3 ) were placed for a length of 30 mm on top of the GBS.
The different experiments are denoted as follows: Exp1 for strong water-wet GBS, Exp2 for weak water-wet GBS, Exp3, Exp4, and Exp5 for the intermediate-wet GBS, and Exp6 for weak air-wet GBS.
A series of 24 fluid invasion experiments were conducted (four replicates for each specific GBS type, Table 2). Air-saturated water was injected with a volumetric flow rate q w = 0.73 mm 3 /s causing a front velocity of 0.007 mm/s. These values correspond to a capillary number of 10 −7 . The water was injected using a high-precision bidirectional syringe pump with a step resolution of 0.046 μm (Fusion 200, Chemyx, Stanford, USA). We note that the accuracy and reproducibility of the experiments depend on this precise flow control at low capillary numbers. We conducted the same experiment using a peristaltic pump and found that typical flow pulsations (pressure fluctuations) during water injection lead to a much higher amount of trapped gas (up to 10%). This is caused by instabilities at the displacement front, which leads to additional bypass trapping.

Computed Microtomography and Image Processing
The fluid-fluid pattern and the pore structure were recorded after each injection experiment and analyzed at three different sections (Figure 6) of the column using computed X-ray microtomography ( Figure 6). An industrial X-ray scanner with a focal spot of 3 μm was used (X-Tek XTH 225, Nikon Metrology, Belgium). The projections were recorded on a Digital PerkinElmer 1620 Flat Panel Detector with a resolution of 2,000 × 2,000 pixels. To maneuver full range of gray values, optimum energy settings (155 μA, 120 kV, and 18.6 W) were used, and beam-hardening artifacts were reduced by using a 0.1-mm copper filter.
The spatial resolution of the CT image (0.013 mm) is important, as it determines the threshold at which a fluid meniscus in a pore throat can be detected. For close random packing, the mean pore-throat diameter can be estimated by d k,min = 0.31 × 0.92 mm = 0.29 mm. Hence, the fluid meniscus can be reasonably resolved by 22 voxels along the diameter.
Image processing was carried out according to the workflow described by Schlueter et al. (2016) using the QuantIm-library (Vogel, 2008) and Fiji ImageJ V.1.50d (Schindelin et al., 2012). The experimental porosity was used as a constraint for selecting appropriate thresholds for the phase segmentation. The raw images were filtered using a nonlocal means filter (Buades et al., 2011). The noise was reduced equally from every dimension (x, y, and z), and contrast at the boundaries between different phases was preserved. Based on the Figure 6. Experimental setup. Sc-1, Sc-2, and Sc-3 show three different sections of the column subjected to μ-CT after each experiment. The column was closed with two-waterproofed caps (a, b) to avoid leakage. A valve was installed at the inlet to control the fluid flow, and a pressure sensor was connected to record water pressure using DASYLab. Blue lines show the connecting PVC tubes.
gray value histogram and threshold detection methods, pore system was separated using the average method (Schlueter et al., 2014). The PSD was analyzed using the maximum inscribed sphere method implemented in ImageJ Plugin BoneJ (Doube et al., 2010). Gas clusters were analyzed by using the labeling method of Hoshen and Kopelman (1976). Segmented images were visualized using VGSTUDIO MAX.

Results and Discussion
In this section, we investigate the percolation transition as a function of the wettability. Therefore, we analyze the front morphology of the invading fluid (Figure 7), the residual gas saturation (Figure 8), the occupation frequency (Figure 9), and the morphology and statistics of the trapped clusters (Figures 10  and 11). Figure 7 shows the percolation transition in the front morphology. The displacement front of Exp1 is compact, due to the strong smoothing effect, whereas the displacement front for Exp5 shows pore-scale fingers and a fractal geometry. The front width is about 2 mm (red dashed window in Figure 7b, bottom figure), which is in agreement with the IP theory (section 2.2.3). Hence, buoyancy has a significant impact on front morphology and leads to a finite correlation length. A representative number of 2-D cuts were used for the calculation of the front width, where boundary effects, for example, wall fingers, were excluded.

Water Resources Research
The displacement fronts of Exp2, Exp3, and Exp4 (not shown in Figure 7) are similar to that of Exp1. The interfacial area of Exp5 is 1.6 times larger compared to that of Exp1. Singh et al. (2017) observed a factor 2 for an oil-water-GBS system. We note that there is some precursor fingering along the Plexiglas wall, because Plexiglas is more wetting compared to the sediment. The optical measurement gives a contact angle of about 70°. This will enlarge the interfacial area but will have only a small impact on bypass trapping inside the sediment. Most of the trapping occurs inside the sediment (Figure 7c).

Capillary Trapping and Occupation Frequency
The smoothing effect has a dramatic effect on the trapping efficiency.
For θ < θ c (Exp1-4), the overlaps at the displacement front, that is, the cooperative front advance, will push away the gas phase (defending phase, Figure 1a). The residual gas saturation is less than 1% for Exp1-4 ( Table 2, first four black points on blue curve in Figure 8).
For θ > θ c (Exp5-6), the front advance is driven by single bursts (Haines jumps) and bypass trapping becomes possible (Figure 1b). This leads to a steep increase in trapping efficiency. The residual gas saturation is about 6% (last black point on blue curve in Figure 8). This steep increase of the residual gas saturation at the critical contact angle θ c = 96°characterizes the percolation transition (blue curve in Figure 8, Cieplak & Robbins, 1990). The trapped gas clusters are shown in Figure 7c.
There is a controversy on whether a second transition from compact displacement to fractal corner flow occurs for contact angles smaller than 45°. Arguments for fractal corner flow are true for 2-D MM experiments (Hu et al., 2018;Zhao et al., 2016), but not for random packs of glass beads. The critical contact angle of 45°is only valid for a rectangle geometry (α = 90°) of a straight flow channel (e.g., a square capillary): θ c = (π − α)/2 = 45°. For a triangular flow channel (α = 60°): θ c = (π − α)/2 = 60°. The pore space in random GBS is (i) interconnected, that is, there is no continuous straight flow channel, and (ii) exhibits an irregular concave-shaped geometry. Hence, the geometrical argument does not hold for GBS. Apart from these simple geometric arguments derived for a straight capillary, the pore space structure (aspect ratio) plays an important role for whether corner flow occurs or not. In Zulfiqar et al. (2020), we studied imbibition into a stochastic quadratic lattice (Si MM, θ ≅ 20°; for details of the MM, see Geistlinger et al., 2019), which should map a straight capillary, perfectly, but we observed no corner flow and no trapping.
In the present study, corner flow was not observed in the high-resolution optical ML experiments (Figures 2 and 3) and in the μ-CT experiments (Figure 7a). The resolution of the μ-CT (13 μm) is sufficient to observe pendular rings (width of the pendular ring w pr ≈ 0.1d grain ≈ 100 μm). If there would be any corner flow, water would be present in pendular rings in front of the main displacement front.
It is instructive to analyze which pores are filled with water and which pores are filled with air/gas. Usually, one expects that air is trapped in the larger pores. However, for θ > θ c air is more wetting than water and should be trapped within the smaller pores. We analyzed the occupation frequency of the pore space and found that indeed air occupies the smaller pores (red curve in Figure 9) and water occupies the larger pores (blue curve in Figure 9). Note that the air occupation frequency is scaled by a factor 8. Figure 8 compares the percolation transition of the gas-water-GBS system (present study) with that of an oil-water-GBS system (red curve in Figure 8, Singh et al., 2017). The oil-water system shows (i) Figure 8. Normalized residual gas saturation S res versus contact angle for gaswater-GBS system (present study; black points and blue fitting curve) and oilwater-GBS system (Singh et al., 2017; gray points and red fitting curve). The best fits to the data are based on an error function fit using the critical contact angle as a fitting parameter. The data are normalized to the maximal residual saturation: 6% for the gas-water system and 15% for the oil-water system. Figure 9. Occupational frequency versus pore diameter for Exp5. The black curve shows the pore size distribution of the total pore space, the red curve shows the pores that are occupied by air/gas (gas cluster size distribution; compare with Figure 7c), and the blue curve shows the water-filled pores. The gas cluster size distribution is scaled by a factor 8.

10.1029/2019WR026826
Water Resources Research a larger critical contact angle (109°), (ii) a weaker slope in the transition region, and (iii) a higher residual gas saturation of 15% compared to the gas-water system. We attribute these differences to the higher viscosity of the oil and the smaller interface tension of the oil-water system. Both physicochemical properties will cause a kinetic slowdown of the front advance. The chance increases that two adjacent menisci will overlap before they burst into new pore throats (Figures 1b and 5a). This causes a higher critical contact angle. Because of the kinetic slowdown effect, the oil phase cannot escape as fast as the gas phase (Figure 1a). This will increase the chance of trapping and could explain the higher residual oil saturation. We believe that the "smooth percolation transition" is also caused by the kinetic slowdown effect. Strictly, this transition is not governed by IP, because kinetic effects are not described by the IP theory (Cieplak & Robbins, 1990).
Calculating the cutoff cluster length caused by buoyancy (equation (12)) for the oil-water system and using a contact angle of 165°( last blue data point in Figure 8), one obtains approximately 5 mm. This length has the same order of magnitude as the sample dimensions (inner diameter of 8 mm) and finite size effects will control the cluster size distribution. Therefore, the experimental cluster size distribution will not exhibit universal scaling, because very large clusters are not present, and hence, the long tail of the power law distribution is truncated by the sample dimension. Maybe this is the reason that the authors could not proof universal scaling for the oil-cluster size distribution. Figure 10 shows the MSE fit (MSE = mean square error estimator) to the volumes of large clusters (last 143 cluster of total 2,742) using the exact scaling exponent τ exact = 2.19 (red curve). The excellent agreement with the universal scaling law predicted by percolation theory (relative data error ε data < 1%) indicates the percolation transition for θ > θ c , that is, IP governs the "local growth modes" (bursts or Haines jumps) at the displacement front, which lead to bypass trapping of the defending phase.

Statistics and Morphology of the Trapped Clusters
Recent studies of cluster size distributions following a power law distribution show that standard MSE fitting is significantly biased compared to maximum likelihood (ML) estimators (Clauset et al., 2009;Goldstein et al., 2004;Iglauer & Wülling, 2016). This bias is often caused by the log-log representation of the data, which gives small clusters (which occur more frequently than large clusters) a high statistical weight. This causes τ values < τ exact and p values < 1 (no linear surface-volume relationship, Figure 6 in . Clauset et al. (2009) developed a statistical framework for discerning and quantifying power law behavior in empirical data. This combines ML methods with goodness-of-fit tests based on the Kolmogorov-Smirnov statistic and likelihood ratios. The corresponding R package "poweRlaw" (Gillespie, 2015) is based on this study. The authors emphasized that the Kolmogorov-Smirnov test is able to detect the best lower limit of cluster size; that is, the data set is reduced as long as the ML objective function shows an absolute minimum. This seems a promising algorithm for selecting the best representative data set. We applied this ML method and obtained a too small lower boundary for the cluster size of 304 voxels (=6 × 10 −4 mm, Cluster Number 400) and a too small τ value of 1.41. This shows that without knowledge about the governing physics, the interpretation of the data is not possible. This is in agreement with the discussion of Stumpf and Porter (2012) about the Critical Truth About Power Laws including the following relevant statements: Figure 10. MSE fits of experimental data (Exp5-sc1, blue points) to a power law distribution (red curve: τ = 2.19). The relative error between data and fitting curve ε data is smaller 1%. Only the large clusters are used for the fit (last 143 of total 2,742 clusters). Figure 11. Linear A i -V i relationship of trapped gas clusters for Exp5.

10.1029/2019WR026826
Water Resources Research 1. Power laws do have an interesting and possibly even important role to play, but one needs to be very cautious when interpreting them. 2. However, a statistically sound power law is no evidence of universality without a concrete underlying theory to support it. Moreover, knowledge of whether or not a distribution is heavy tailed is far more important than whether it can be fit using a power law.
Based on the experimental evidence that the percolation transition occurs for θ > θ c , we search for the physical lower cluster size at which the reduced data set shows universal scaling. If we apply the discrete ML estimator (poweRlaw software) to the large clusters data set shown in Figure 10, we obtain the universal scaling exponent τ = 2.13 (relative error ε τ = 3%).
The second characteristic of the IP theory is a prediction about cluster morphology. The surface of large clusters is linearly correlated to their volume (section 2.2.2). The data for Exp5 (θ = 99°) shown in Figure 11 exhibit a strong linear A i -V i relationship with a regression coefficient R 2 = 0.99. This regression coefficient near 1 is obtained for all column sections and all four replicates demonstrating again the percolation transition that occurs at θ c = 96°.

Conclusions
Our main conclusions are as follows: 1. For the gas-water-GBS system the wettability is a crucial parameter that governs the percolation transition at θ c = 96°. A small increase in θ (about 10°) causes a qualitative change in the dynamic pattern formation. Such a jump-like, highly nonlinear behavior is typical for a first-order thermodynamic (static) phase transition. This is the reason that Cieplak and Robbins (1988) called this transition a dynamic phase transition. This transition occurs within 10°compared to the broad transition region of about 100°of an oil-water-GBS system. We attribute this to a kinetic slowdown effect, due to the higher viscosity of the defending fluid oil (2 orders of magnitude) and its lower surface tension compared to the gaswater-GBS system. 2. This crossover in pattern formation causes a crossover in trapping efficiency; that is, it changes from no trapping to maximal trapping. 3. The fluid occupation frequency of pores clearly demonstrates that the gas clusters are trapped within the small pores. This trapping behavior is typical for the more wetting fluid. 4. For θ > θ c , our experimental results show that IP governs the fluid displacement. Strong indicators are the universal scaling behavior of the size distribution of large clusters (relative data error ε data < 1%) and their linear surface-volume relationship (R 2 = 0.99). We note that the large clusters exhibit a cutoff length of about 2 mm caused by buoyancy that is much smaller than the column diameter of 19 mm. We think that this is the reason that the occurrence of large trapped gas clusters is statistically representative.
The key question is whether this interesting dynamical phase transition will occur in (i) sediments with rough surfaces, for example, sands and sandstones, and (ii) in sediments with more irregularly shaped grains. Preliminary results in Si MMs with irregular grains, a smooth surface, and a contact angle of about 40°show residual gas trapping of about 60%. Previous studies on capillary trapping in natural sands (estimated contact angle about 20°;  show residual gas trapping of about 15%. Future research will investigate the impact of grains with rough surface and irregular geometry on the Cieplak-Robbins transition and will be published in Part 2 of the study (Zulfiqar et al., 2020).