Synchrotron-based pore-network modeling of two-phase flow in Nubian Sandstone and implications for capillary trapping of carbon dioxide

Depleted oil fields in the Gulf of Suez (Egypt) can serve as geothermal reservoirs for power generation using a CO2-Plume Geothermal (CPG) system, while geologically sequestering CO2. This entails the injection of a substantial amount of CO2 into the highly permeable brine-saturated Nubian Sandstone. Numerical models of two-phase flow processes are indispensable for predicting the CO2-plume migration at a representative geological scale. Such models require reliable constitutive relationships, including relative permeability and capillary pressure curves. In this study, quasi-static pore-network modelling has been used to simulate the equilibrium positions of fluid-fluid interfaces, and thus determine the capillary pressure and relative permeability curves. Three-dimensional images with a voxel size of 0.65 micro m3 of a Nubian Sandstone rock sample have been obtained using Synchrotron Radiation X-ray Tomographic Microscopy. From the images, topological properties of pores/throats were constructed. Using a pore-network model, we performed a sequential primary drainage-main imbibition cycle of quasi-static invasion in order to quantify (1) the CO2 and brine relative permeability curves, (2) the effect of initial wetting-phase saturation (i.e. the saturation at the point of reversal from drainage to imbibition) on the residual-trapping potential, and (3) study the relative permeability-saturation hysteresis. The results improve our understanding of the potential magnitude of capillary trapping in Nubian Sandstone, essential for future field-scale simulations. Further, an initial basin-scale assessment of CO2 storage capacity, which incorporates capillary trapping, yields a range of 14-49 GtCO2 in Nubian Sandstone, Gulf of Suez Basin.

Schematic diagram of key processes of CO2 sequestration in a saline aquifer. As caprocks are typically not perfectly horizontal, the injected CO2 will tend to move updip and once the injection of CO2 ends, the entire CO2 plume will typically move updip. At the trailing edge of the CO2 plume, this leads to re-imbibition of native formation fluid (e.g., brine), reducing the CO2 saturation in the pore space. When this imbibition process is completed, some amount of CO2 remains trapped by capillary forces. This "residual trapping" effectively immobilizes a large portion of the injected CO2 and can, under favorable geologic conditions, remove the necessity of a sealing unit (caprock) to ensure permanent CO2 storage underground. Depths are not to scale.
[B] Typical pore-scale configurations of fluid distributions at various pore-space saturations.

Pore-scale trapping mechanisms
Numerical modeling of two-phase flow in a natural porous medium at the pore scale is often employed to obtain two-phase constitutive relations such as capillary pressure and relative permeability curves, and to understand pore-scale displacement mechanisms (e.g. Patzek, 2001;Blunt et al., 2013;Berg et al., 2016, and references therein). The application of high-flux synchrotron radiation in the hard X-ray range provides improved time and spatial resolution, permitting the investigation of the internal geometry and topology of microstructures down to sub-micrometer levels (Donoghue et al., 2006). The constant improvement in both imaging capacity and computational power allows usage of discretized images of the pore structure as a computational mesh to perform fluid-flow simulations directly (Blunt et al., 2013;Gostick, 2017). Pore-scale models can be broadly divided into direct numerical simulation models (such as volume-offluid and Lattice Boltzmann methods), and pore-network models (PNM). In addition to being computationally cheap, PNM can simulate flow in pore networks that preserve the pore connectivity expressed in pore-space images. Each (network) element is assigned a set of geometric properties (e.g. size and shape) that closely match the values measured directly from the images.
Using a high-performance computing cluster, the multiform geometrical and essential topological properties of the Nubian Sandstone pore network (PN) have been directly extracted from Synchrotron Radiation X-ray Tomographic Microscopy (SRXTM) data. Further, we have performed quasi-static porenetwork modeling, based on the intricate structure of the Nubian Sandstone PN. In this work, we are interested in the flow processes occurring between the two immiscible fluid phases at the trailing edge of the CO 2 plume and these processes' relevance for trapping CO 2 in Nubian Sandstone.
To investigate imbibition processes at the trailing edge of the CO 2 plume, the following pore-scale processes are considered: 1. Piston-like (pl ) advance is characterized by the displacement of the nonwetting phase from a narrow pore throat into an adjoining pore body that was initially filled with the wetting phase. It causes a drop in capillary pressure, but it does not result in the trapping of the nonwetting phase.
2. Snap-off (so) occurs when two wetting layers in a pore body swell, such that they touch and coalesce, disconnecting the non-wetting phase within the neighboring pore bodies from the main connected pathway.
3. Co-operative pore-body filling (cpf ) refers to the invasion of a non-wetting-phase-filled pore throat by arc menisci (AM) or wetting layers initially present in corners, crevices and rough surfaces of the pore body. At a critical capillary pressure, the AMs fuse and the center of the throat spontaneously fills with the wetting phase. This can lead to the trapping of the non-wetting phase.

Objectives
The main objective of this paper is to characterize the two-phase properties of the Nubian Sandstone as it pertains to CCS and CPG by implementing a quasi-static PN model. Of particular interest is the quantification of the potential for immobilization of CO 2 by capillary forces and how this compares to other sandstones. In addition, the PN model is used to simulate sequential primary drainage/main imbibition cycles to obtain capillary pressure-saturation and relative permeability-saturation curves. Such constitutive relationships are important inputs for large-scale numerical simulations of CCS/CPG operations. These results complement rigorous laboratory experiments for understanding fluid-displacement mechanisms in porous media.
We have endeavored to replicate pressure, temperature, and salinity conditions that are similar to those found in the reservoir of interest. An important outcome of this study is the basin-scale assessment of CO 2 storage capacity for Nubian Sandstone in the Gulf of Suez basin in which capillary trapping is considered.

Rock formation characterization
The Nubian Sandstone -a major reservoir-rock formation is found in Egypt and extends regionally-has the potential to serve as a deep geological formation for long-term CO 2 storage and related applications, such as CPG-based power generation, as discussed in the introduction.  presented a three-dimensional static reservoir model for Nubian Sandstone, which includes the formation's geology and rock physics properties. For the current study, Nubian Sandstone blocks have been retrieved from the central Gulf of Suez (cGOS).
We conduct a routine optical microscopy analysis, employing an automated QEMSCAN (Quantitative Evaluation of Minerals by Scanning Electron Microscopy) Quanta 650F technique at the University of Geneva, Switzerland, providing quantitative mineralogical compositions at a scanning resolution of 5 µm. Here, the mineral identification maps are based on back-scattered electron values, energy-dispersive X-ray spectra, and X-ray count rates. These tools are integrated to assess the mineralogical composition and rock type.
The QEMSCAN analysis shows that the majority of the distinguished mineral phases in the Nubian Sandstone are medium to fine quartz grains (85%) with some kaolinite (8%), 5% k-feldspar, and 2% other minerals, trapped within the pore system, as shown in Figure 2A. The quartz grain types (which, in turn, are mainly monocrystalline with sharp edges; Figure 2B) are cemented by silica. This Palaeozoic quartz arenite sandstone was deposited in continental to fluvial braided river systems. It represents the primary prolific conventional reservoir in the cGOS. Such a quartz-dominated system suggests that neither mineral dissolution/precipitation nor clay swelling require consideration during the quasi-static simulations.
For the laboratory rock physics experiments, three cylindrical plugs of 25.5 mm diameter and > 30 mm height were cored from the rock block parallel to each other to investigate the intrinsic rock properties and quantify possible anisotropies. A lathe was used to produce a high degree of parallelism between the upper and lower surfaces of the cylindrical plugs, with ±5 µm precision. The plugs were dried in an oven at 100 • C for at least 48 hours to remove trapped pore fluid. Before performing further measurements, the plugs were stored in a desiccator with silica gel to prevent moisture from re-entering. Non-peer reviewed preprint: April 16, 2020 Using a gas-displacement-based apparatus (He-pycnometer Accupyc II 1340, micromeritics), the matrix volumes of the cylindrical plugs were determined under ambient conditions in a range that keeps the standard deviation of the measurements below 5%. Beforehand, the plugs' masses were measured, using a high-precision balance (±5.0 ×10 −7 kg tolerance). The bulk volume, matrix volume, dry mass, grain density, and effective porosity (i.e. the volume of interconnected pores) were calculated. A 2D slice through the SRXTM-volume raw dataset of the Nubian Sandstone (voxel size is 0.65 µm 3 ). Dark gray areas are pore space (air), while light gray areas represent the mineral grains (quartz). Due to the experiment setup, the data that are used for pore-network modeling are located only within the marked white square. The corresponding file name of this slice 'SRXTM_Nubian_Sandstone_06_rec.8bit_0084.tif' is given in the full raw dataset. The full SRXTM raw dataset (2560×2560×4320 pixels and 8-bits) is provided as a supplementary to this publication.
A cylindrical subplug (with 6 mm diameter and ∼ 10 mm length, that was recovered from the main plug) was used to perform high-pressure Mercury Intrusion Porosimetry (MIP) to evaluate the pore radius, pore-size distribution (PSD), and porosity-related characteristics of Nubian Sandstone. The PSD observations are based on the behavior of a non-wetting liquid in a capillary, following the Young-Laplace relationship. The experiment was set up initially at low pressure (40 MPa), followed by a stepwise increase in pressure up to 400 MPa. At 400 MPa, mercury can be forced to invade pores as small as 5 nm in radius, constituting microporosity (Kärger, 2011). At high capillary pressures, the wetting phase (air) is pushed farther into the pore corners, as the radius of curvature of the interface decreases.
To validate the results of our PN quasi-static simulator, we compare the pore-volume density data and the scaled air-mercury P c −S w experiment curve with the simulation results extracted from the 3D rock model. Mercury Intrusion Porosimetry (MIP) shows a lognormal pore-size distribution of the Nubian Sandstone samples. MIP also shows that the majority of the pore bodies fall within the resolvable range of our imaging capacity ( Figure 3).

Imaging technique
High-resolution Synchrotron Radiation X-ray Tomographic Microscopy (SRXTM) can be employed to investigate in a non-destructive and quantitative way the internal topology of micro-structures. The SRXTM measurements were conducted on a cylindrical subplug (with 6 mm diameter and ∼ 10 mm length that was recovered from the main plug). The parallelism of the ends of the cylindrically shaped subplug was achieved by gently grinding the cuttings, first on the side on a rock saw blade, then by hand, using sandpaper (grit 120).
We performed the imaging work at the TOmographic Microscopy and Coherent rAdiology experi-menTs (TOMCAT) beamline station, at the third-generation synchrotron facilities, at the Swiss Light Source, Paul Scherrer Institute, Switzerland (see Marone and Stampanoni, 2012, for technical specifications). For the best-contrast phases (pore and solid), the sample was exposed to a parallel beam of monochromatic synchrotron X-ray radiation for 500 ms exposure time at 26 keV photon energy and 401.25 mA ring current. This beam energy is optimized to provide a sufficient photon flux, capable of penetrating the 6 mm diameter subplug and ensuring the best image contrast (signal-to-noise ratio). The lower-energy X-rays that hit the cylindrical subplug and do not improve imaging were filtered out using 100 µm Al, 10 µm Cu, and 10 µm Fe sheets, while the remaining X-rays hit the cylindrical subplug. The field of view (FOV) covered 1.3 × 1.4 mm 2 of the original size of the subplug.
The transmitted X-rays were converted into visible light by a cerium-doped lutetium aluminium garnet (LuAG : 20 µm) scintillator and projected at 10× magnification onto a high-speed CMOS camera (PCO.Edge 5.5; PCO AG, Germany) with 2560×2560 pixels, leading to an effective pixel width of 0.65 µm. A sample-to-scintillator distance of 36 mm yielded a small amount of edge enhancement in the images.
Each tomogram was computed from 1801 projections (500 ms exposure time) over a 180 • rotation by a gridded Fourier transform-based reconstruction algorithm with a Parzen filter using a filtered backprojection algorithm (Marone and Stampanoni, 2012). Projections were magnified by microscope optics and digitized by a high-resolution CCD camera, which results in 8-bit (256 gray values) tiff-format images. We provide, through the online access of this article, the SRXTM images data set along with a summary of the experiment conditions and the characteristic properties of the SRXTM images.
The reconstruction center was found for the first and last image in the sequence and linearly interpolated between these two values for the others. The reconstructed volumes were filtered with a 3×3×3 median filter, segmented with local connectivity-based thresholding, and processed further as described in the next section. The image datasets were then cropped around the plug so each one consisted of around 2160×2160×4320 ∼ 2.02 ×10 10 pixels with a voxel size of 0.65 µm 3 at 10-fold optical magnification. A 3D volume with dimension 2160×2160×4320 voxels was cropped from the original 3D raw data set to avoid poor image quality around the edges of the scan volume.

Image processing and segmentation
The gray-level SRXTM images of the Nubian Sandstone have only one major identifiable mineral phase (Quartz; from Section 2) with sharp edges ( Figure 2B). The attenuation coefficient of the X-ray (grayscale intensities) in the mineral phase is higher than the attenuation coefficient in pore spaces (darker gray). A bimodal histogram for each image, therefore, was applied using Otsu's method of automatic thresholding (Otsu, 1979); Figure 4A. The threshold value was set for the local minimum to alleviate spurious creation of catchment basins.   Figure 2B by implementing the multi-level Otsu method (Otsu, 1979). [B] Application of the local Thickness filter with a spherical structuring element to get actual information about the pore-size distribution (i.e. determination of the radius of the largest sphere that fits inside the pore space; Rmax) from the CT images.
[C] The pore-size distributions determined using the simulated Local Thickness filter. The pore-volume distribution indicates that approximately 80% of the total pore volume comprises pores with a radius smaller than 55 µm. The histogram of the pore distribution shows that the majority of the pore radii are in the range of 17 − 50 µm, closely matching the experiment values obtained from Mercury Intrusion Porosimetry; see Section 2 and Image segmentation is done by delineating and labeling sets of pixels within the domain of graylevel SRXTM images into a number of disjoint non-empty sets of objects depending on the gray-level histogram, neighborhood information, and prior knowledge of the image statistics (Soille, 2014). The watershed transformation, which treats pixel values as local topography, was used to segment the images. However, watershed segmentation does not straightforwardly provide acceptable results due to oversegmentation. Alleviating this effect will be elaborated upon in the following. The Local Thickness filter replaces every voxel in the binary images with the radius of the largest spherical kernel that fits inside the void space ( Figure 4B). Analysis of the histogram of the voxel values provides information about the pore-size distribution as illustrated in Figure 4C.
After thresholding, a labeling procedure was performed on the obtained binary images to isolate and quantify pore-space voxels from solid voxels using the Euclidian distance map. The distance transform of a bi-level image is defined by the shortest distance from every (pore-space) pixel to the nearest non-zero foreground-valued (solid) pixel and is employed to distinguish between pore bodies and pore throats. Every voxel is assigned to a pore, such that the distance map increases in the direction of the pore centre; ( Figure 4D). The algorithm also finds throat surfaces by looking for restrictions in the void space that bound pores: here the distance map increases on either side of the throat surface. The distance map is inverted as the watershed grows towards high intensity.
A Gaussian blur filter was applied to the distance transform map to smooth it and to minimize spurious peaks caused by the flat nature of the solids voxelated in the images. In order to minimize over-segmentation, a manual optimization yielded 0.35 and 5 as optimal values for the Gaussian blur filter and the Local Thickness filter, respectively.
Removal of the erroneous peaks that lead to over-segmentation was done individually and in an iterative way, using a minimal cubic structuring element during morphology dilation. The similarity between the maximized dilated peaks and the smooth distance map allows the recognition of new peaks. The comparison between the old peak(s) and the new peak(s) enables eliminating spurious errors. This is followed by applying the marker-based watershed transformation at the local minimum of a hypothetical 3D topography, from which basins are flooded to segment the images and separate the pores ( Figure 4F). It is possible to compute properties of pore-network elements, including diameters, coordinates, surface area, and volume, for pore bodies as well as diameters, perimeters, and lengths for pore throats. A detailed description of the extraction algorithm and the parameter definition of the PN elements is given in Gostick (2017).

Pore-network construction
Since the pore morphology of Nubian Sandstone is quite complex ( Figures 2B and 4A), its topological quantification is challenging. Its 3D intricate pore network was constructed using the discretized elements, resulting from the marker-based watershed segmentation method (see Section 3.2 and Figure 4). In this paper, the geometrical properties of the extracted pore network refer to the distributions of the pore-body and pore-throat radii, lengths, and volumes. Additionally, relevant topological traits include the spatial connectivity, which prescribes how pore bodies are connected via pore throats (i.e. coordination number, z) and a shape-describing parameter for each pore-network element that will be used during numerical fluid-flow modeling ( Figure 5). Coordination numbers range from 1−39 for the whole extracted network, with a mean coordination number of z ≈ 4.1, which is close to the value found by Dong and Blunt (2009) for Fontainebleau and Berea sandstones. Isolated clusters of pores (i.e. pores with a coordination number of zero) were removed from the pore network.
A three-dimensional visualization of the extracted pore network is presented in Figure 5A, where the red spheres represent pore bodies and gray cylinders represent the pore throats, where both appear to be properly located within the rendered voxel image. The regularity of the pore-element shapes is of course an over-simplified representation of the natural porous medium, lacking information of high relevance for two-phase flow. Many studies have focused on idealizing the geometry of pore-network elements to mimic such structure complexities and angularities (Mason and Morrow, 1991;Patzek and Kristensen, 2001;Øren and Bakke, 2003). These idealized pore elements will be discussed in more detail in Section 4.1. The distributions of inscribed radii for the extracted pore-network elements are shown in Figure 5, while Table 1 provides the statistical information of the pore-network properties, used in the simulations.

Representative elementary volume
Before simulating two-phase fluid flow within the constructed pore network (PN), the existence of a representative elementary volume (REV) has to be confirmed by the spatial distribution of its components Non-peer reviewed preprint: April 16, 2020   . The voxel SRXTM image shows the solid phase in gray, while the extracted pore elements are overlain. The pore bodies and pore throats are respectively rendered as red spheres and gray cylinders that are smaller than their actual sizes to improve visualization.
[B] Histograms of pore-body and pore-throat radii as well as pore-throat lengths for the extracted pore network of the Nubian Sandstone.
[C] The distribution of the coordination number, (z), in the Nubian Sandstone pore network. Table 1: Parameter statistics of the pore-network structure for the Nubian Sandstone sample, extracted by the stochastic pore-network extractor, and used in quasi-static PN modeling. Network size (domain) mm 3 is 1.41×1.41×2.82.

Pore-network parameters
Pore bodies Pore throats a Extended diameter of spherical pore element is defined as the maximum value of the global distance map lying within each pore region. b Inscribed diameter of spherical pore element is obtained in the same manner as the extended diameter, but a local distance map of just the pore region is used. c Cross-sectional fraction was calculated after determination of the dimensionless shape factor (G) (derived from Mason and Morrow, 1991); see Section 4.1 for more details. The shape-factor concept is used to describe the geometry of the pore-network elements and to relate it to the displacement processes being simulated.
Non-peer reviewed preprint: April 16, 2020 (e.g. solid-void spaces) to satisfy a lower and an upper bound of the fluid-flow field. These bounds are required to ensure that the REV is sufficiently larger than a single pore and smaller than the flow domain, such that the phase saturations are independent of the REV size and representative of the macroscopic domain. Figure 6A shows that the physical properties (here the porosity) are highly fluctuating as the sub-element volume considered is small, and the values stabilize and approach the experiment values, measured via different techniques and for different sizes, as the averaging volume increases. The minimum size of the REV for Nubian Sandstone was found to be ∼ 10 9 µm 3 ( Figure 6A). In addition, Figure 6B shows variations in the porosity profile as a function of distance across two opposite sides of the sample. These variations are very small and are thus neglected. The two sides correspond to the inlet and outlet boundaries in our PN simulations (Section 4.2). Moreover, it is important to note that the anisotropy parameter of permeability, i.e. the ratio between the vertical and the horizontal permeability, is less than 5%. Therefore, its significance, in relation to fluid flow effects can be neglected during our quasi-static simulations. We employ direct numerical simulations (DNS), based on the Lattice Boltzmann method (LBM), to determine the permeability anisotropy in the three orthogonal directions (Appendix A.2). In summary, we consider here that the Nubian Sandstone sample is mostly homogeneous in nature and that therefore the local fluid flow vectors always point in the opposite direction of the potential gradient during the simulations. The DNS results show no preferential flow direction, which, judging from the homogeneous nature of the samples (by visual inspection), is not surprising.

Numerical methods
Over the past few decades, several studies have estimated fluid-flow properties from digital images of rocks (e.g. Valvatne and Blunt, 2004;Blunt et al., 2013;Krevor et al., 2015;Rasmusson et al., 2016Rasmusson et al., , 2018Qin and van Brummelen, 2019). In this paper, the laws of mass and volume conservation (incompressible and immiscible two-phase, laminar flow in a natural porous medium) are considered, where one fluid, supercritical CO 2 (scCO 2 ), is the nonwetting phase, while the other (brine with a salinity of ∼ 40 % ) is the wetting phase. The thermophysical properties of these fluids are listed in Table 2. The salinity value refers to the surface water salinity in the Gulf of Suez after Gavish (1974). The assumption here is that seawater was injected into the reservoir in the past to enhance the recovery of oil. We assume that initially, the brine fully filled the pore network of the Nubian Sandstone.
A quasi-static PNM simulator for scCO 2 invasion after Qin et al. (2016) has been modified for this study, yielding capillary pressure and relative permeability curves: • The capillary pressure curve is achieved by simulating an invasion-percolation process.
Non-peer reviewed preprint: April 16, 2020 • Two-phase relative permeabilities are calculated at each saturation state of the two-phase quasistatic displacement process, once system equilibrium is reached.
The parameterization of the capillary pressure curve is done with the well-known van Genuchten (VG) model (Equation 1), given as (van Genuchten, 1980): In Equation (1), S e is the effective saturation of the wetting phase, where S w , S wr , and S nwr denote the wetting phase saturation attained by the primary drainage curve, the residual saturation of the wetting phase, and the residual saturation of the nonwetting phase, respectively, and α d,i Pa −1 , n, and m = 1 − 1/n are the VG fitting parameters (Table 3). The subscripts 'd' and 'i' refer, respectively, to the drainage and imbibition process.  (2010)). c Density was obtained after Span and Wagner (1996)'s equation of state.

Shape factor
A 2D slice of Nubian Sandstone ( Figure 2B) shows highly irregular geometries of the pore elements. Even though it is not feasible to reproduce these geometries explicitly in the PNM, it is important to consider their effect on the displacement process. It is necessary, therefore, to replace this geometry by an equivalent idealized geometry which leads to mathematically tractable problems, while at the same time reproducing the main features of the pore space, pertinent to two-phase flow. Mason and Morrow (1991) introduced the shape-factor concept to describe the geometry of porenetwork elements (i.e. pore bodies and pore throats) and their angularity for reliable two-phase flow simulations, where the wetting phase occupies the angular corners and the nonwetting phase represents the bulk (Figure 7). The most common approach is to choose the shape -i.e. either a circle, a triangle, or a square (C−T−S) -whose shape factor (as defined in Equation 2) is closest to that of the real pore, as the idealized pore geometry, where A represents the cross-sectional area of a pore or throat and P is its corresponding perimeter. Our marker-based watershed segmentation provided: (1) the region volume, (V ), by summing the number of voxels in that region; (2) the extended radius, (R ext ), given as half of the maximum value of the global distance map, lying within each pore region, which may lead to an overlap of two pore bodies; and (3) the inscribed radius, (R ins ), given as half of the maximum value of the local distance map (i.e. Euclidean distance for a single pore space), lying within each pore region. The length of an idealized pore element, (L ipe ), is the summation of the extended radius and the inscribed radius (i.e. Equation 3). The cross-sectional area of an idealized pore element, (A ipe ), and its perimeter, (P ipe ), can, therefore, be defined by Equations (3) and (4), respectively, given by: Non-peer reviewed preprint: April 16, 2020 For pore throats, the marker-based watershed segmentation counts the number of voxels in the throat with a distance transform value of 1, then multiplies the result by the voxel length to provide the pore throat perimeter. Its inscribed radius can be found from the maximum of the global distance transform. With Equations (3) and (4), and substitution with Equation (2), we obtain a new equation that describes the shape factor of three-dimensional idealized pore elements, From Equation (5), the shape factor values of √ 3 36, 1 4π, and 1 16 correspond to the 2D crosssections of a scalene triangle, square, and circle, respectively. Statistical properties of pore-network elements using the newly reformulated shape factor (Equation 5) are given in Table 1. The details of how to calculate the corner half-angles β 1 , β 2 , and β 3 , with the convention of [0 ≤ β 1 ≤ β 2 ≤ β 3 ≤ π 2], as illustrated in Figure 7, for a given shape factor, are given in Patzek (2001).
Figure 7: Wetting-phase (gray) distribution in a pore element, which has a triangular cross-section. The three half-corner angles, β1, β2, and β3, are 15 • , 30 • , and 45 • , respectively. θa is the static advancing contact angle, R is the radius of the inscribed circle, r is the radius of the nonwetting-wetting-interface curvature, and b denotes the distance from each apex to the corresponding nonwetting-wetting interface along the solid wall (red).

Quasi-static pore-network modelling
The pore network (PN), initially fully filled with the wetting phase, is brought in contact with the nonwetting and wetting reservoirs at the inlet and outlet boundaries, respectively. The other boundaries are set to wall (no-flow) boundary conditions. Both reservoir pressure values are set to zero initially. When the primary drainage starts, we gradually increase the nonwetting reservoir pressure by small increments (less than 5 Pa). After each increase, we check which pore bodies are invaded by the nonwetting phase. A pore body will be invaded by the nonwetting phase, once the following two conditions are both satisfied: • The nonwetting reservoir pressure is larger than its entry pressure, • Both nonwetting and wetting phases in the pore body have pathways to the inlet and outlet reservoirs, respectively.
Then, the phase occupation in both pore bodies and pore throats is updated under the equilibrium condition. The capillary-pressure relationships, with respect to the wetting saturation used, are given in Appendix A.1. The primary drainage process is terminated once a pre-set maximum drainage pressure is achieved. At this stage, the wetting phase is present in the tight pores and crevices of the network, due to pore-structure trapping, and resides in the corners of the surface roughness of the pores. We use an interfacial tension of σ = 0.0202 N/m. For main imbibition, we decrease the nonwetting reservoir pressure step-wise, by a small deduction (less than 5 Pa), allowing the wetting phase gradually to fill the pore bodies in each step. After each decrease, we check the scenarios of cooperative pore-body filling (cpf ) and snap-off (so) in pore bodies and pore throats. For the cooperative filling of a pore body, if the local capillary pressure is larger than the threshold capillary necessary to fill the pore body, and the nonwetting and wetting phases have pathways to the inlet and outlet reservoirs, the pore body will be fully filled with the wetting phase by spontaneous imbibition. When snap-off happens in a pore throat, the pore throat will be completely occupied by the wetting phase. When snap-off happens in a pore body, the nonwetting phase will be locally trapped. For the remaining pore bodies and pore throats, the phase occupation is updated by the piston-type movement of arc menisci (AMs) or the main terminal meniscus (MTM) filling (Ma et al., 1996) under the equilibrium condition. The relationships for checking snap-off events and the threshold capillary pressure for the cooperative filling, used here, are given in Appendix A.1.
Finally, in both drainage and imbibition processes, as long as the phase occupation needs to be updated, a search operation needs to be conducted to check if the nonwetting and wetting phases exhibit pathways that are connected to the inlet and outlet reservoirs, respectively. If both pathways can be found, the phase occupation will be updated to the new one. Otherwise, the update is aborted due to trapping.

Calculation of relative permeability
For a given stage in the drainage or imbibition process, the conducting pathways of the nonwetting and wetting phases are known. For either phase, its volumetric conservation in a pore body is given by: Here, i is the pore-body index, ij is the pore-throat index, N i is the coordination number, κ is the phase indicator, P [Pa] is the pressure, and K κ ij m 3 /Pa/s is the conductivity, calculated by where l [m] is the length, and g m 4 /Pa/s is the conductance, given in Appendix A.1. With proper inlet and outlet pressure boundary conditions, Equation (6) can be used to solve for the pressure field. Then, Darcy's law is used to obtain the relative permeability for the given network saturation value.

Results
In the following sections, the pore-space geometry of the Nubian Sandstone sample, as determined experimentally using the relevant data from the MIP experiment, serves as an input for the image-based PNM extraction approach. In all the simulations that follow, we run primary drainage-main imbibition cycles using the quasi-static pore-network model, while varying the saturation at the point of reversal from drainage to imbibition.

Pore-network calibration
We assume that the pore throats are volumeless elements of the pore network and account for resistance to flow, while pore bodies account for the pore volume but do not provide any resistance to flow. On a local level, however, the extracted PN of the Nubian Sandstone contains slate-shaped pores that have an aspect ratio (i.e. ratio of pore-body to pore-throat sizes) smaller than one. Only small spheres can be inscribed inside pore spaces with such pore-to-throat aspect ratios, but the throats could still be quite large. In this case, such a division into pore bodies and throats has little physical meaning. Therefore, we adjust the pore-throat radii until the measured capillary pressure-saturation relationship, obtained from the MIP experiment, is fulfilled. This is done by assuming that the aspect ratio of any pore-body/porethroat pair must be greater than, or equal to, a given threshold value. The resulting minimum aspect ratio, R asp , for the analyzed Nubian Sandstone sample in this work is 2.27, compared to 2.25 for Berea Sandstone (Øren and Bakke, 2003). This calibration is done manually by trial-and-error. The P c − S w relationship used for the calibration is obtained by extending the results of the MIP experiments from the mercury-air fluid pair, using a Leverett-J scaling (Leverett, 1941) to the CO 2brine system (Table 3). In addition, a van Genuchten (VG) model was fitted to the P c − S w curve as shown in Figure 8. (1) Measured via He pycnometry on a cylindrical plug of the same rock (25.4 mm diameter cylindrical sample). This value is based on averaging three plugs.
(4) The permeability was estimated, based on steady-state, single-phase flow across a constant pressure difference between the inlet and the outlet of the pore network.  Table 3 summarizes the intrinsic hydraulic and capillary properties of the Nubian Sandstone, using experimental and computational methods. The results of the quasi-static single-phase flow simulation shows that the absolute permeability is in agreement with the value obtained from the MIP experiment. It can be seen that the absolute permeability of the network 2.00 ×10 −12 m 2 is close to that of the plug sample MIP experiment; 1.70 ×10 −12 m 2 . The permeability obtained from direct numerical simulation, using the Lattice Boltzmann method (our in-house code, LBHydra, Saar (2010a,b, 2013)) is with 2.57 ×10 −12 m 2 slightly higher.

Residual trapping in Nubian Sandstone
In this section, we investigate pore-scale fluid-displacement and trapping mechanisms by simulating primary drainage-main imbibition cycles. The simulations use the calibrated pore network of the Nubian Sandstone sample (Section 5.1). We assume water-wet conditions and that the contact angle is constant for each process, i.e. that the receding contact angle, θ r , is 0 • (primary drainage), and the advancing contact angle, θ a , is 35 • (main imbibition). The primary drainage simulations stop when the pore-space saturation of the non-wetting phase reaches a target value, S nw,i . This is the initial condition for the subsequent main imbibition process, where the wetting saturation increases till the non-wetting phase is no longer connected to its reservoir. A visualization of the primary drainage-main imbibition cycle, with S nw,i = 10% and an advancing contact angle of 35 • , is shown in Figure 9. The brine-filled pores are shown in blue, while red represents CO 2 -filled pores. The figure shows clearly that most residual trapping takes place under low capillary pressures, i.e. towards the end of the imbibition cycle. In the simulations, S nw,i is varied between 0.1 and 0.9, in order to assess the effect of initial saturation on the residual nonwetting saturation. Further discussions on the relationship between S nw,i and S nw,r (i.e. the initial-residual relationship, IR) are provided in Section 5.4. Furthermore, we investigate two distinct processes, whose competition influences the residual trapping at the pore scale. The first is changing the system wettability from strongly water-wet (i.e. θ a = 0 • ) to more neutrally water-wet (i.e. θ a = 50 • ). The second is the sensitivity of residual trapping to the minimum aspect ratio, R asp . In addition to R asp = 2.27, obtained during the pore-network calibration, we chose R asp values of 1.14 and 4.54 to analyze statistically the fluid distribution, the saturation map, and the trapping events. Figure 10 shows, how the amount of residually trapped CO 2 differs with contact angle and minimum aspect ratio, R asp . Generally, more CO 2 trapping occurs at lower contact angles and higher minimum aspect ratios. However, this does not hold true for the lowest minimum aspect ratio of R asp = 1.14. In this case, S nw,r increases with increasing contact angle, and even surpasses those of the higher aspect ratios at high contact angles.  Figure 9: Snapshots of the quasi-static simulations, with a wetting-phase (brine) saturation, at the point of reversal from drainage to imbibition, of Snw,i = 0.1. The brine-filled pores are shown in blue, while red represents CO2-filled pores. The upper panel shows a representative process of primary drainage at three different capillary pressures, while the lower panel shows imbibition. We use only half of the extracted pore network, shown in Figure 5A, as it is more appropriate for the simulations to have a cube-shaped domain.  , as a function of the minimum aspect ratio of the pore network. The figure shows that an increase in the amount of trapped CO2 corresponds to a decrease in the contact angle for the medium (2.27) and high (4.54) minimum aspect ratios. In contrast, the amount of residual CO2 trapping increases with increasing advancing contact angle for the smallest minimum aspect ratio (1.14).

Relative permeability
In this study, the predicted values for the relative fluid permeabilities in the Nubian Sandstone sample are obtained via quasi-static pore-network modeling. These values are then compared to the two-phase experiment results, reported in Reynolds et al. (2018) for the Bunter Sandstone (UK) as well as to the VG functions, which were fitted to the capillary pressure curve, obtained from our MIP experiment ( Figure 11). Generally, there is a good agreement between the simulated relative permeabilities and the calculated values from the VG function (fitted with the capillary pressure curve). The Nubian Sandstone relative permeabilities tend to be higher than those observed experimentally for the Bunter Sandstone, but the same trend is seen. The end-point nonwetting-phase relative permeability, at the irreducible water saturation for the Nubian Sandstone sample, obtained from our PN simulations, is k r,nw (S w,ir = 0.099) = 0.794.
:  We investigate the dependence of hysteresis, of the relative permeability curves, on the pore-scale phenomena described in Section 5.2 (i.e. wettability and minimum aspect ratio), as shown in Figure A.4. Based on the observations from our pore-network study, we infer that the hysteresis in the k r,nw − S w relationship is primarily caused by the dependency of the swelling of arc menisci in pore corners during imbibition. The lower the initial wetting saturation (i.e. the wetting saturation when drainage switches to imbibition), the more the nonwetting phase (here CO 2 ) in the pore bodies tends to get disconnected at the end of the imbibition process, thereby becoming more isolated and immobilized, in part via cooperative pore-filling.

Initial-residual saturation relationship for Nubian Sandstone
Two-phase modeling at the continuum scale can benefit significantly from realistic residual trapping (IR) constitutive relationships in order to quantify the magnitude of trapping (Krevor et al., 2015). Several empirical trapping models have been developed to describe residual trapping. The most widely used model for water-wet systems is that of Land (1968). In this model, the residual saturation is calculated by where, C is the Land trapping coefficient, which is found by fitting a curve to a plot of the residual saturation of the nonwetting phase (S nw,r ) versus the initial saturation at the end of drainage and prior to imbibition (S nw,i ). The higher C is, the lower is the amount of trapped nonwetting fluid (here CO 2 ), and therefore the weaker the hysteresis between drainage and imbibition. We use the calibrated pore network (Section 5.1) to estimate the trapping curve for the CO 2 − brine pair in Nubian Sandstone and compare it with the trapping curves that were experimentally determined for similar rock types by Akbarabadi and Piri (2013) and Reynolds et al. (2018).   Figure 12 shows that the capillary trapping curves increase monotonically with initial saturation. At S nw,i ∼ 0.6, the relationship is nearly linear. Interestingly, Nubian Sandstone has relatively good capillary trapping characteristics with up to 45% of the CO 2 initially present remaining as a trapped fluid after imbibition. A similar residual trapping value has been observed by Spiteri et al. (2008) for Berea Sandstone. It is worth noting, however, that in a mixed-wet system, residual trapping would be greatly reduced (Spiteri et al., 2008).

Discussion
Overall, the Nubian Sandstone sample is relatively homogeneous, with only small variations in the porosity profile along its length ( Figure 6). The developed quasi-static pore network (PN) simulator demonstrates the sensitivity of residual trapping in such a homogeneous sandstone to the key parameters wettability and minimum pore-body-to-throat aspect ratio.

Sensitivity analysis of Pore-Network modelling
Displacement mechanisms: Generally, three different pore-scale displacement mechanisms are possible during imbibition, as explained in Section 1.2. Piston-type displacement is always favored over cooperative pore-body filling and snap-off, wherever piston-type displacement is topologically possible. However, in order for nonwetting-phase trapping to occur, cooperative pore-body filling or snap-off is required. Figures 13A and 13B show two histograms of the number of occurrences of each displacement mechanism in the pore bodies and pore throats, respectively. One can see that, for the two higher minimum aspect ratios, co-operative pore-body filling (cpf ) is the main mechanism responsible for the trapping of the nonwetting phase in the pore bodies. cpf occurs most often for R asp = 2.27. Snap-off (so) is all but absent in the pore bodies. For R asp = 1.14, however, so is slightly more common than cpf . In addition, the number of snap-off events in the pore throats increases with increasing minimum aspect ratio. Such behavior is expected, as it is easier for the wetting phase, present in the corners of the pore throats, to connect (leading to snap-off) in narrow throats. Effect of wettability on residual CO 2 trapping: For R asp = 2.27 and R asp = 4.54, the relationship between the amount of residually trapped CO 2 and the contact angle is straightforward. The lower the contact angle, the higher the entry pressure required for a nonwetting CO 2 bubble in a pore body to invade a (narrow) pore throat. Hence, as shown in Figure 10, S nw,r reduces with increasing contact angle. This behaviour is consistent with previous experiment results, using core-flooding computerized tomography on different sandstones (Chatzis et al., 1983;Herring et al., 2013;Geistlinger et al., 2014;Krevor et al., 2015).
However, this tendency is inverted for R asp = 1.14. It has been shown in imbibition experiments with low-aspect-ratio capillary micromodels, that strongly water-wet displacement yields better displacement efficiencies (sweep) than mixed-wet systems (Mattax and Kyte, 1961). Hence, lower contact angles result in smaller residual trapping values in this case. The reason for this effect is that systems with high capillarity tend to show fluid distributions that are more homogeneous (Mattax and Kyte, 1961), while larger clusters of nonwetting fluid can persist when capillary forces are less relevant (Figure 14), yielding more residual trapping.

Non-peer reviewed preprint: April 16, 2020
Effect of the pore-body-to-throat aspect ratio on residual trapping: For strongly water-wet systems, higher minimum aspect ratios yield more trapping (Figure 10). This is again caused by the higher entry pressures that narrower pore throats have. As has been pointed out throughout this discussion, the case of R asp = 1.14 behaves decidedly differently than the other two cases with larger aspect ratios. It should be noted though that such a low aspect ratio is uncommon in natural rocks. As such, the trapping behavior and trends depicted by the cases with R asp = 2.27 and R asp = 4.54 are more representative for CO 2 residual trapping in sandstone.

Implications for CO 2 storage in the Gulf of Suez Basin
We explore the implications of this research, as a first attempt, for future implementation of CCS in the Egyptian context. Equation (9) is used to estimate the geologic CO 2 storage capacity. It is consistent with methods used in previous studies to assess the prospective geologic storage of buoyant fluids in subsurface formations (van der Meer, 1995;Doughty et al., 2001;Kopp et al., 2009;Goodman et al., 2011;NETL, 2015), where, M eff CO2 is the effective storage capacity [kg], V bulk res is the bulk reservoir volume m 3 , and ρ CO2 is the CO 2 density kg.m −3 as a function of the corresponding reservoir temperature and pressure. The dimensionless CO 2 -storage efficiency factor, ζ eff [−], represents the fraction of the total pore volume that can be occupied by the injected CO 2 . Doughty et al. (2001) proposed a method in which ζ eff is estimated by a combination of the parameters provided in Equation (10).
Here, φ avg is the average formation porosity [−], C g and C h are coefficients for the geometric capacity [−] and the heterogeneity capacity [−], respectively. The coefficient C i is intrinsic capacity [−], which is controlled by multi-phase flow and transport phenomena; Equation (11).
where C ig and C il are the intrinsic capacity [−] for the CO 2 and brine phase, respectively. C ig can be estimated using Equation (12).
where, the S g is the average saturation of CO 2 within the CO 2 -plume. We use the residual CO 2 saturation value of 0.4 obtained from the current pore-network analysis. Meanwhile, the C il coefficient can be estimated using Equation (13).
where, S l = 1 − S g is the average brine saturation [−], χ g l is the average CO 2 mass fraction dissolved in the brine phase [kg/kg], and ρ l /ρ g stands for the density ratio between the CO 2 and brine phases [−]. We employ the commonly used thermodynamic models to describe the mole fraction (solubility) of CO 2 for a CO 2 -brine system provided by Duan and Sun (2003) and Duan et al. (2006). The densities, as a function of the relevant pressure and temperature, are numerically obtained using the Span and Wagner (1996) equation of state for CO 2 and Driesner (2007) for brine (considering the brine salinity discussed above).
In the Gulf of Suez basin, Nubian Sandstone is subdivided into three main zones with high reservoir qualities and interfingering of two relatively thin-shale aquitards between the sandstone aquifer zones. Such layer-type heterogeneities tend to enhance sequestration capacity by counteracting gravitational forces (Doughty et al., 2001). Therefore, we choose a heterogeneity capacity coefficient C h of 0.5 (van der Meer, 1995). A range of geometrical capacity coefficients C g is taken from Kopp et al. (2009): 0.19-0.63. Stratigraphically, Nubian Sandstone overlies the basement complex throughout most of the Gulf of Suez basin with an average thickness of 0.45 km (Alsharhan and Salah, 1997). The depth map for the basement complex was interpreted from the aeromagnetic and the geological data by Mesheref et al. (1976);Farhoud (2009). Therefore, the depth to the top of the Nubian Sandstone sequence was constructed by adding a vertical offset of 0.45 km to the interpreted depth map of the basement complex. The resultant map is shown in Figure 15A with a maximum depth of 4.67 km.
Temperature and pressure generally increase with burial depth but have opposite effects on CO 2 density. Following the geothermal gradient of 35.7 • C/km and pressure gradient of 10.18 MPa/km after  [B] The distribution of in-situ CO2 densities across the Gulf of Suez. CO2 density is calculated at the equivalent temperature and pressure after Span and Wagner (1996)'s equation of state. The red outline bounds the potential area (with 9 950 km 2 footprint) in the Gulf of Suez for CO2 to be in the supercritical phase.

Parameter
Value Model et al. (2006) ρg kg/m 3 647.7 a Span and Wagner (1996) ρ l kg/m 3 959.47 a Driesner ( a These values are calculated as a function of brine salinity, and reservoir temperature and pressure (Hefny, 2020). b The area estimation (marked with the red outline in Figure 15B) bounds the potential area in the Gulf of Suez where the depth to Nubian sandstone will be > 1.5 km and the CO2 will be in the supercritical phase. The depth map is constructed based on the interpretation of aeromagnetic and geological data by Mesheref et al. (1976) for the basement rocks and modified after Farhoud (2009). c Volume calculation is based on an average thickness of 0.45 km after Alsharhan and Salah (1997).
Non-peer reviewed preprint: April 16, 2020 , the range of CO 2 density was found to be 640 -624 kg/m 3 . Hefny (2020) estimated χ g l [−] and ρ l /ρ g ratio [−] as 0.057 and 1.48, respectively, both are calculated for an average temperature (118.24 • C) and pressure (34.7 MPa), and brine salinity (0.66 mol/kg; Table 2) found in the Gulf of Suez. All calculated coefficients and parameters used to estimate the CO 2 storage capacity in Nubian Sandstone are provided in Table 4. In conclusion, our ζ eff estimation for Nubian Sandstone was found to be in the range of 0.005 -0.017 [−].
The M eff CO2 estimation covers a 9 950 km 2 surface area for the Gulf of Suez (marked by the red outline in Figure 15B), where the subsurface conditions allow CO 2 to be in the supercritical phase. We conclude that the Nubian Sandstone in the Gulf of Suez can store up to 14 -49 GtCO 2 , which alone can meet the Egyptian needs to store 250 MtCO 2 emitted per annum (Crippa et al., 2019) for the upcoming >227 years.

Conclusions
• We investigate geologic CO 2 storage in a highly permeable formation, typical for Nubian Sandstone, Gulf of Suez (Egypt). This investigation consists of (1) constructing a realistic 3D pore network model that represents the characteristic features of Nubian Sandstone, (2) developing a quasi-static pore-network numerical simulator that mimics in-situ conditions that are similar to those prevailing at a CO 2 -storage site, at the trailing edge of the CO 2 plume, and (3) predicting the two-phase flow characteristics.
• Two-phase constitutive relationships, including the capillary pressure and relative permeability curves, are computed for water-wet rocks at a low capillary number. We determine a Land trapping coefficient of C = 1.2 for the Nubian sandstone sample. These relationships can be employed in field-scale numerical models to estimate the extent of CO 2 -plume migration, at a representative geological scale.
• The estimated capillary-trapping curves for Nubian Sandstone are in good agreement with those observed experimentally for similar rocks. This confirms, that residual trapping due to hysteresis can be a key mechanism for long-term CO 2 storage in Nubian Sandstone.
• One aspect, which has not been taken into account in this study, is contact-angle hysteresis, due to the lack of experiment measurements of the contact angle in Nubian Sandstone. We, therefore, consider uniform wettability in this work. Wettability hysteresis may be a topic of interest for future studies.
• Our sensitivity study of the impact of wettability (i.e. advancing contact angle) on the magnitude of residual trapping shows that the trapped amount of CO 2 increases with decreasing contact angle for the cases, where the minimum pore-body-to-throat aspect ratio is 2.27 or 4.54. The inverse is the case when the minimum aspect ratio is only 1.14.
• The pore-network model developed in this work improves our understanding of the different trapping mechanisms in Nubian Sandstone, as they pertain to CCS-and CPG-related applications.
• The result of a conservative approach for a basin-scale estimation of CO 2 storage capacity shows an attractive potential to store 14 -49 GtCO 2 in Nubian Sandstone, Gulf of Suez basin. This assessment is based on a horizontal caprock and incorporates the arrangement of several injection wells accessing all fault-blocks with a compartment-like basin. The existence of many of these required wells, due to the locally active hydrocarbon industry, is another attractive factor, contributing to potentially making CCS and CPG techno-economically feasibile in the Gulf of Suez.

Supplemental material
Supplemental material is provided to enable the reproducibility of our results. In addition, the SRXTM raw dataset of the Nubian Sandstone sample (in tiff format) has been made publicly accessible for the scientific community through the ETH Zurich Research Collection repository [DOI: https://doi.org/10.3929/ethzb-000377881], which also contains descriptions of imaging conditions and characterizations of the rock samples through laboratory measurements .
Non-peer reviewed preprint: April 16, 2020 There are two imbibition mechanisms for a pore body, namely, filling from arc menisci (AMs) and the main terminal meniscus (MTM). For the MTM filling, the capillary pressure (i.e., the entry pressure) will keep constant as the increase of wetting saturation until a pre-set value of saturation is reached.
In this work, we set it to 0.96. Then, a transitional function of the wetting saturation is used to force the capillary pressure to fast approach zero. Figure A.1 shows three types of functions. The type 1 is primarily used in this work. The capillary entry pressure for a pore body or a pore throat is calculated by the following nonlinear equations based on the MS-P theory (Patzek, 2001): where, A eff is the effective cross-sectional area occupied by the nonwetting phase, A is the cross-sectional area of the pore body, L ns and L nw denote the length of nonwetting-wetting interface and the length of the nonwetting-wall interface, respectively, and for the remaining notations, one can refer to Figure 7.
For a pore body with the cross section of the right triangle described in Figure 7, assuming R = 20 [µm] and the surface tension (σ) of 0.073 [N/m], the obtained entry capillary pressure versus the static contact angle is plotted in Figure A.2A. Also, the corresponding distances b 1 , b 2 , and b 3 are plotted. The distances go to zero if corner flow cannot be formed.
For the imbibition filling from AMs, taking the parameters used in Figure A.2A, the capillary pressure versus saturation curves are plotted in Figure A.2B for four contact angle values. Notice that each curve is truncated at its snap-off saturation. The snap-off event is triggered once more than one AM touch each other. For a triangle cross-section, we check the following three inequalities: where I denotes the side length. If one or more are false, the snap-off is triggered. When the snapoff occurs, the nonwetting phase is trapped in the pore body. As mentioned in Section 5.1, no volume is assigned to pore throats. The capillary pressure in a pore throat is approximated by the capillary pressure in its connected upstream pore body.

A.1.2 Phase conductance
Relations of capillary pressure and phase conductivities for idealized pore elements are shown in Figure  7. Flow resistance in both pore bodies and pore throats is taken into account. If a pore is fully filled by either nonwetting or wetting phase, the single-phase conductance is given by Patzek and Kristensen (2001): where, g is the phase conductance.
If a pore is occupied by both nonwetting and wetting phases, the conductance for the wetting phase is calculated as follows. First, for an AM (see Figure 7, with b = 1, the dimensionless corner wetting area is given as proposed by Patzek and Kristensen (2001): Further, we define the shape factor for the corner wetting as: Then, the semi-empirical dimensionless conductance with the perfect slip boundary condition between wetting and nonwetting phases is calculated by Patzek and Kristensen (2001): 1 4 π −G w + 2 ln S w (A.12) Finally, the dimensional one is obtained asg w = b 4 µ g w . Taking the parameters used in Figure A.2. Figure A.3 shows the curves of conductance versus saturation for four contact angle values. Obviously, for the same saturation value, a larger contact angle gives larger wetting-phase conductance along the corners. For the nonwetting-phase conductance, it is approximated by Qin (2015):

A.1.3 Relative Permeability hysteresis
The hysteresis patterns appear in relative permeability curves are shown in Figure A.4. These simulated k r − S w curves are presented as a function of initial wetting-phase saturation, wettability and aspect ratios.

A.2 Direct Numerical Simulation (DNS); Absolute Permeability
At the pore-scale, the permeability and flow field has been determined using Lattice Boltzmann Method (LBM) to solve the conservation equations in the complex geometry of Nubian Sandstone. A representative geometry of the Nubian Sandstone was reconstructed from a high-resolution Synchrotron Radiation Xray Tomographic Microscopy (SRXTM) binarized volume of 1.4×1.4×1.3 mm 3 . For the stability of LBM simulation, we performed a voxel-averaging in order to down-scaling the resolution SRXTM's images to 270×270×270 pixel while keeping the 3D internal geometrical structures of unchanged. We used ImageJ open-source image analysis tool to lowering the voxel size to the value of 5.13 µm 3 . See the supplementary section for accessing the raw SRXTM grey images (original pixel size). Fluid flow inside the reconstructed Nubian Sandstone is performed simulating a single phase flows using an in-house LBM numerical permeameter (LBHydra; Saar, 2010a,b, 2013) Non-peer reviewed preprint: April 16, 2020 the estimated permeability of pore-network modeling. The D3Q19 model is implemented in the solution of 3D incompressible fluid flows in the LBM. At the solid boundary of the pore space, no-slip and no-flow BCs can typically be utilized by "bounce back" implementation. During the LBM simulations, once flow was deemed to have reached steady state ( Figure B.1A), the permeability of the medium was calculated based on the Darcy's law and after the conversion from lattice units into physical units be achieved. Figure B.1B displays the 3D distribution of streamlines of velocity in z-direction of Nubian Sandstone. The pressure difference between the inlet and outlet, which is expressed by the density difference in the simulation, is 0.0005 in lattice units. As shown in Figure B.1, the streamlines represents the continuous pathways from the inlet to the outlet are quite straight, resulting in quite high permeability of Sandstone. The permeability is determined in the three flow directions at the same scanned data of Nubian Sandstone sample to investigate its permeability's anisotropy ( Figure B.1B and Table B.1). Based on the coupled LBM method, the predicted permeability results are 2.38 × 10 −12 , 1.62 × 10 −12 , and 3.71 × 10 −12 that corresponding, respectively, to permeability m 2 in x-, y-, and z-directions. The mean value of the permeability is 2.57 × 10 −12 m 2 . The standard deviation is 1.06 × 10 −12 m 2 . Its anisotropy, the ratio of the horizontal to the vertical permeability, to be < 1 confirming the homogeneity of Nubian sandstone.  Figure 1A to show the geographic locations of the measured samples around Gulf of Suez (Egypt). The orientation of the drilled plugs and its abbreviations is visualized in Figure 1B. Henceforth and due to the high homogeneity of the collected Nubian Sandstone rock sample, the reset of the dataset will include the experiments and simulation only for NS06 sample (grey shadow in Table 1, the closest outcrop sample to the studied reservoir in Hefny et al., 2020).

Porosity directory
Folder name ! "Porosity_and_Permeability" Subfolder name ! "Helium_Porosity" Files names ! "Sample-name_orientation.txt" File 02 name ! "porosity_density_calculations.csv" Table 1: Geographic locations and their coordination for the collected Nubian sandstone rock samples around Gulf of Suez ( Figure 1). In the table below, the gray-shaded row contains the sample "NS06" which has been used to perform the current work.  [right] Orientation of the cored plugs relative to bedding of the main oriented rock block. The 25.5 mm plugs show the drilling directions of the directional measuring of rock physics properties. Consequently, we used the following convention: NS stands for Nubian sandstone abbreviation following by X, Y, I or Z and corresponding to the measurements in the four mutual directions. The subplug with 6.5 mm diameter is the sample which is used for synchrotron CT scanning.

Permeability measurements
This folder contains the permeability measurements following transient-step method up and downstream reservoir pressure curves vs. times only for one Nubian Sandstone. The gas-permeability measurements for Nubian Sandstone plug (NS06 sample; three different directions) were performed using an oil-based hydrostatic pressure medium apparatus developed in-house, following the transient pulse decay method (Brace et al., 1968). Each file tabulates up-and down-stream reservoir pressure values vs. experimental times.

Mercury Intrusion Porosimetery
This folder represents the Mercury Intrusion Porosimetry (MIP) measurement, which is performed at ClayLab, Institute for Geotechnical Engineering, ETH Zürich, Switzerland. The MIP measurement records the capillary pressure incremental vs. the volume of the non-wetting phase increase. The experimental conditions are also provided.

Petrographic analysis
This folder contains QEMSCAN (Quantitative Evaluation Materials Scanning Electron Microscopy) surface area mineralogy for Nubian Sandstone, which is performed at Department of Earth Sciences, Geneva University at a scanning resolution down to 5 µm.

Images acquisition at the Swiss Light Source Synchrotron
The ultrahigh-resolution 3D volumes of Synchrotron radiation X-ray tomographic microscopy (SRXTM) are provided in this folder (in .tif format), together with descriptions of imaging conditions; Table 2.
The SRXTM dataset has the original (uncompressed) gray images of Nubian Sandstone acquired at the TOMCAT beamline of the Swiss Light Source Paul Scherrer Institute, Switzerland, with resolutions down to 0.65 µm. The size of the image dataset exceeds 26 GB (total of 4320 .tif images). In order to avoid uploading-downloading difficulties (such as browser and internet connection), we split it into 6 subunits and each subunit contains 720 .tif images.

TOMCAT imaging directory
File 01 name ! "TOMCAT_Experimental_Conditions_NS06_Sandstone.txt" Folders names ! "Synchrotron_CT_images_part_number" File 01 name ! "SRXTM_info_range_number.txt" Images names ! "SRXTM_Nubian_Sandstone_06_rec.8bit_number.tif" Each excel file contains one sheet for a complete cycle of primary drainage at Sw,i of 0.1 as well as at least 7 sheets for main imbibition under different Sw,i of 0.1 -0.9.

Reservoir Properties for Nubian Sandstone
A three-dimensional compartmentalized reservoir model is undergoing by the authors for Nubian Sandstone in the Gulf of Suez (Egypt). Simplified reservoir properties are provided in Table 3.