Towards Large-Scale DRP Simulations: Generation of Large Super-Resolution images and Extraction of Large Pore Network Models

Representativity and accuracy of digital rock physics (DRP) simulations depend strongly on the size of the image volume and the resolution obtained. Even with one of the fastest DRP simulation techniques like pore network modelling, simulation volumes have typically been limited to few cubic millimetres for highly resolved images. In this paper, a super-resolution technique named enhanced super-resolution generative adversarial network (ESRGAN) is used to obtain well-resolved images with large field of view and to generate micro-CT images with resolution enhancement factors of × 4 and × 8. Subsets of resulting ESRGAN images were tested against the same volume (of acquisitions resolved at high and low resolution) by comparing petrophysical properties of interest. Pore network extraction and multiphase simulation results showed that pore size distribution, porosity, permeability, drainage capillary pressure and relative permeability curves obtained using ESRGAN images were more accurate. Large images, however, pose subsequent limitations on DRP simulations as pore network extraction code needs a lot of memory to process them (usually more than 60 GB of RAM for 15003 voxels image). Thus, we present a novel stitching strategy that is developed to enable the extraction of pore networks on such large images. Several validation cases of this method are presented to test the accuracy of the results from stitched networks on single- and multiphase flow properties. Finally, our stitching tool was used to generate two large networks of 3.6 million and 9.2 million elements, respectively, from two large ESRGAN images of approximately 49003 voxels. Application of ESRGAN on micro-CT images with resolution improvement factors × 4 and × 8. Testing the image quality enhancement by performing single- and multiphase flow simulations. Generation of large super-resolved images with size up to 4688 × 5160 × 4800 voxels. Development of a stitching methodology to extract pore networks from these large images. Application of ESRGAN on micro-CT images with resolution improvement factors × 4 and × 8. Testing the image quality enhancement by performing single- and multiphase flow simulations. Generation of large super-resolved images with size up to 4688 × 5160 × 4800 voxels. Development of a stitching methodology to extract pore networks from these large images.


Introduction
Achieving large-scale digital rock physic (DRP) simulations is a challenging subject. A trade-off is often needed between the field of view and the resolution of the micro-CT images. A small field of view questions the representativity of the computed properties, while a larger field of view (at low resolution) impacts the precision of simulation results (Saxena et al. 2017). To circumvent this compromise, an evolving strategy is to digitally enhance the resolution of 3D micro-CT images using a deep learning approach so as to obtain micro-CT images with large field of view and good resolution (da Wang et al. 2020;Chen et al. 2020;Ahuja et al 2022). In this work, we use a deep learning approach called enhanced super-resolution generative adversarial network (ESRGAN).
Recent advances in deep learning methods have led to major advances in super-resolution imaging. The first family of super-resolution methods is based on convolutional neural networks and is named SRCNN (Dong et al. 2015;Lim et al. 2017;da Wang et al. 2019;Yu et al. 2018). They have been used by several researchers to improve photographic and medical images. Recently, these methods have been applied to micro-CT rock images and have proved very successful in capturing the edges of the pores as they compute a pixelby-pixel loss. However, SRCNN tends to decrease the pixel-by-pixel difference during the training and increase the peak signal-to-noise ratio (PSNR). Unfortunately, higher PSNR does not necessarily reflect the perceptually better results (Ledig et al. 2017). Ledig et al. (2017) proposed a super-resolution generative adversarial network (SRGAN) for which they used a deep residual network with skip-connections and defined a new perceptual loss using high-level feature maps of the VGG 19 network (Simonyan and Zisserman 2014). This resulted in a better capture of the texture of the images. da Wang et al. (2020) applied SRGAN to micro-CT images of rocks and obtained better visual check with the high-resolution images as the texture was recovered. Chen et al. (2020) proposed a cycle-consistent generative adversarial network CycleGAN-based super-resolution approach that they applied on real rock images and showed that their method improved the porosity estimation on the considered acquisitions. Wang et al. (2018) improved the SRGAN and introduced an upgraded method that they named enhanced super-resolution generative adversarial network (ESRGAN). In fact, they changed the basic building unit of the generator network and introduced residual-in-residual dense block (RRDB) without batch normalization as they observed that these batch normalization layers introduced some undesirable artefacts. They also used relativistic average GAN (Jolicoeur-Martineau 2018) instead of the classical GAN allowing the discriminator to predict the relative realness instead of the absolute value. Finally, the perceptual loss has been improved and was computed using the feature maps before activation, and this provided better brightness consistency and texture recovery. Overall ESRGAN provided better results than SRGAN and won the first place in the PIRM2018-SR challenge (Blau et al. 2018). It is worth noting that most existing super-resolution studies mentioned above work with synthetic images. But recently, Ahuja et al. (2022), have developed a new super-resolution method called Siamese-SR model that focused on edge-semantics. They applied this method to improve the resolution and texture of micro-CT rock images and compared their results to other existing methods by applying SRGAN, ESRGAN, EDSR and SPSR methods to the same rock images. While their simulation results for Siamese-SR based porosity, permeability and MICP were very encouraging and outperformed other methods, the study was only demonstrated for a small resolution enhancement factor of 2 (in each direction). Furthermore, the study did not assess the performance of these methods on relative permeability curves, which are one of the most desired petrophysical parameters to be predicted by DRP simulations. In this work, we have demonstrated the application of ESRGAN method on micro-CT rock images for larger enhancement factors of 4 and 8 and used these results to evaluate the performance of ESRGAN in predictions of porosity, permeability, capillary pressure and relative permeability curves with respect to the original high-resolution (HR) and low-resolution (LR) counterparts.
Several simulation techniques have been developed to study single-and two-phase flow in porous media. Direct numerical simulation (DNS) methods such as finite volume method and lattice Boltzmann have been widely used to compute petrophysical properties of rocks directly on micro-CT images (Alpak et al. 2018;Raeini et al. 2012;Yang 2013). However, these methods are computationally expensive for two-phase and capillary dominated flow regime, which characterize most oil reservoirs. Pore network modelling (PNM) technique, on the other hand, simplifies the pore structure and idealizes its geometry through a pore network extraction (Dong and Blunt 2009) that is later used in a pore network simulator to compute single-and two-phase flow properties. The latter technique is very efficient computationally and can be used to simulate larger volumes (Regaieg and Moncorgé 2017). However, the older models (Oren et al. 1998;Valvatne and Blunt 2004) used geometry simplifications to compute the node-to-node conductivities, which is a critical parameter in petrophysical properties computations. In order to take the best of both approaches, Raeini et al. (2017) performed a single-phase DNS simulation on the rock image and coupled it with a pore network extraction. This results in a pore network with conductivities computed for the original geometry.
Although pore network modelling is computationally efficient, for the very large images obtained from super-resolution, memory swiftly becomes a limiting factor. Indeed, the extraction code needs large amount of memory to extract networks from large images (usually more than 60 GB of RAM for 1500 3 voxels image). To perform simulations in large volumes, Jackson et al. (2022) acquired low-resolution images, applied super-resolution to several sub-volumes in which they performed pore network simulations and finally upscaled the flow using continuum-scale simulation. However, pore network simulations in each grid could be biased by finite size effects if they are not large enough (Regaieg et al. 2021). Therefore, extraction of large pore networks is important even when the simulations are coupled with an upscaling strategy. Several researchers have developed strategies to overcome the extraction limitation. Rabbani et al. (2019) divided the image into several non-overlapping subdomains and extracted a pore network from each subdomain. These networks were then stitched together to get a pore network on the full image. The stitching was done by analysing the interaction of the throats at the intersection of two neighbouring subdomains. Although their algorithm improved the computational efficiency, the final pore network was found to differ from the network extracted from a single large domain. In fact, the rules used to connect the different subnetworks were introducing uncertainty. Similarly, Kohanpur and Valocchi (2020) developed a pore network stitching approach where pores were generated from statistics on the boundaries between the subdomains in order to connect them. Again, this can introduce uncertainties and inconsistencies between the reference pore network and the stitched pore network. Khan et al (2020) parallelized the watershed step and then performed network extraction on the entire image. This helped them to have a gain of 50% of memory usage and to extract a pore network from 2000 3 voxels image. However, they used overlap to deal with interface between the subdomains, and evaluating the needed overlap was difficult to manage. This methodology does not solve completely the memory issue as the extraction step is still performed on the full image and is not parallel. In this work, we have developed a new stitching methodology where the different extractions are performed on overlapping subdomains, thus removing the need to generate pore elements statistically to connect the different domains. Furthermore, we have decoupled the extractions and we were able to extract a pore network from 4688 × 5160 × 4800 voxels image.
This paper is structured as follows: First the ESRGAN method implemented in this work, our training strategy and our newly developed stitching method are described. Subsequently, quality of generated super-resolution images (generated with image enhancement factors of 4 and 8) is tested with respect to their high-resolution (HR) and low-resolution (LR) volume counterparts by assessing petrophysical metrics obtained from single-and multiphase flow simulations. Following validation of ESRGAN results, the developed stitching methodology is tested for sensitivity and consistency. Finally, two large-scale super-resolved images (up to 4900 3 voxels) are generated, and respective pore networks extracted using our stitching strategy. Primary drainage simulation results for these large networks (up to 9.2 million total elements) are shown.

ESRGAN: Method and Training Strategy
In this work, we have implemented enhanced super-resolution generative adversarial network (ESRGAN) proposed by Wang et al. (2018). This method is based on two phases of training: (1) PSNR step and (2) GAN step. In the PSNR step, L1 loss function is minimized during the training of the generator. In this stage, the borders of the pores are improved, but the texture of the rock is not captured. The weights of the trained generator are used as the starting point of the second stage of the training that we call GAN training. Pre-training with PSNR helps the GAN to have more visually pleasant results. The loss function of the second stage has three terms: the adversarial loss, the perceptual loss and the content loss. For the adversarial loss, Wang et al. (2018) proposed to use a relativistic average GAN (Jolicoeur-Martineau 2018) where loss function does not optimize discriminator to distinguish data real or fake. Instead, RaGAN's discriminator distinguishes that "real data is not like average fake data" or "fake data is not like average real data". According to Wang et al. (2018), this helps to learn sharper edges and more detailed textures. The perceptual loss helps the network to capture perceptually relevant differences like textures or the clays in the rock images. These are determined from features maps obtained from VGG19 before the activation. We do not use the weights of the pre-trained VGG19 network, instead we re-train it with rock micro-CT images in order to have more relevant features and we have observed an improvement of the results after that. Finally, the content loss is the L1 pixel-by-pixel between the generated and the high-resolution images.
The loss for the generator is therefore: where L percep is the perceptual loss, L 1 is the 1-norm distance between recovered image G(xi) and the ground truth y, and λ (equal to 5 ×10 −3 ) , η (equal to 1 ×10 −2 ) are the coefficients to balance different loss terms. We use the network architecture introduced by Wang et al. (2018) with several residual dense residual in blocks (RDDB) and without batch normalization layers. We use the variant with 23 RDDB as it helps achieving better trainings. We note that we have changed the up-sampling layer and we introduced sub-voxel convolution instead of the deconvolutional layers in order to reduce checkerboard artefacts as suggested by Yu et al. (2018). We have prepared two versions of our neural networks one with resolution improvement ×4 and another one with improvement ×8 . We highlight that we do not aim to train a general neural network that can be applied to different rocks after a single training. In this work, our objective is to use a neural network on the same rock it was trained on in order to maximize our chances of having more accurate results. Having this in mind, it is important to have a fast training in order to be able to handle new rocks when needed. Therefore, distributed data parallel (DDP) in Pytorch was used to parallelize the training with multi-GPUs multinodes parallelization and we have used 36 GPUs in parallel in the trainings presented in this paper (taking 2 days per training). Multiprocessing with distributed data parallel duplicates the model across multiple GPUs, each of which is controlled by one process. During training, each process loads its own minibatches (6 images in each) from disk and passes them to its GPU. Each GPU does its own forward pass, and then, the gradients are allreduced across the GPUs. Gradients for each layer do not depend on previous layers, so the gradient all-reduce is calculated concurrently with the backwards pass to further alleviate the networking bottleneck. At the end of the backwards pass, every node has the averaged gradients, ensuring that the model weights stay synchronized. The generator is first trained by minimizing the L1 loss; we call this step the PSNR training and it allows to recover the edges of the pores but not the texture as we can see in the example in Fig. 1. The second stage of the training is the GAN training where more complicated features like the clays and the textures are recovered as we can see in the example in Fig. 1.
We have observed that the larger the crop size the better the performance of the training. Therefore, we used 384 × 384 voxels for the high resolution and 96 × 96 in ×4 variant and 48 ×48 in the ×8 variant. We used two acquisitions of the same volume with different resolutions (LR and HR) for both our trainings and for each training 3000 crops are prepared as our training data.

Pore Network Stitching
Pore network models do not involve two-phase direct flow simulation in a 3D CT image or a reconstructed digital rock. Instead, it needs a pore network that is extracted from 3D reconstructions. Various algorithms for network extraction exist to extract the skeleton of the 3D model that carries the essential geometric and topological information of the underlying pore system. In this work, we use a pore network extraction platform called GNextract (Raeini et al. 2017). GNextract is first used to reconstruct an upscaled version of the 3D segmented image of a rock in the form of a network of pore elements where the single-phase flow conductances in each pore are derived by solving the Stokes equation in the original geometry using OpenFOAM. Unfortunately, the extraction code needs large amount of memory to extract large images (more than 60 GB of RAM for 1500 3 voxels image). Furthermore, although OpenFOAM simulations are parallelized, they are computationally expensive and require a large number of resources and/or simulation time for large images. Therefore, to overcome these limitations, a novel stitching process has been developed on networks extracted from overlapped sub-volume of a given image. This process, unlike pore network extraction, does not take a voxel image as an input, but only pore networks. Figure 2 illustrates the stitching process and shows that the stitched network is formed by a first part conserved from PNM1, a second part conserved from PNM2 and a third part obtained on a transition zone where rules inspired from pore network extraction method in general are defined to choose the elements to conserve-some elements being obtained from PNM1 whilst others from PNM2. It is important to have a transition zone as there could be some inconsistencies induced by boundary effects between both PNMs for the same part of the overlapping sub-image.
The key part of pore network stitching is to define rules to choose between elements from PNM1 or PNM2 in the transition zone. To define these rules, elements of both PNM1 and PNM2 within the transition zone are gathered into families of pore bodies and pore throats. On the one hand, pore bodies are defined as belonging to the same family if their embedded spheres intersect each other. On the other hand, families of pore throats are defined by gathering all throats connecting the same families of pore-bodies. In the end, Fig. 2 An illustration of the pore network stitching process only the largest element in each family is conserved in the stitched network (as illustrated in Fig. 3).
If a large extraction is required, several stitching operations are performed, thus reducing the memory usage and accelerating the computations.

Factor 4 Enhancement
We tested the method on two rocks: a Bentheimer and a Reservoir Sandstone that we call Reservoir C (Fig. 21). For each rock, the data were split into training, cross-validation, and test set. For the Bentheimer rock, images of two samples (S1 and S2) were available and both samples were scanned fully at 5.8 microns resolutions and a smaller volume of each was scanned at 1.45 microns. Therefore, the first sample (S1) images were divided into training and cross validation sets and the second sample (S2) images were used as a test set. We observed that ESRGAN was sensitive to the histogram of the images. Therefore, the histogram of S2 LR (test) images had to be matched to the histogram of S1 LR (training) before applying the generator.
For the inference, an up-sampling with a bicubic interpolation in the Z direction is performed on the low-resolution dataset to reduce the pixel size in the vertical direction. Then, the generator is applied slice by slice on 2D images with 796 × 820 pixels to generate 3D image of 3184 × 3280 × 12168 with a voxel size of 1.45 microns. We noticed that as the Fig. 3 Illustration of stitching rules for choosing pore bodies (up) and pore throats (down) in the transition zone, where red elements belong to PNM1, while blue ones are related to PNM2 low-resolution training images were also upsampled in the Z direction, the neural network learns to correct for the blurriness created from the up-sampling step. After that, a test subset from the super-resolution image is considered and compared to the high-resolution zoom of the same volume.
For Reservoir Rock C, similar methodology was used. The main difference was having a single sample scanned fully at 5.36 microns resolution and a zoom of a sub-volume was performed with 1.34 microns. Again, we divided the dataset into training, cross-validation, and test datasets. After the training and the resampling of the low-resolution dataset in the Z direction, the generator is applied on 2D slices of 1172 × 1290 pixels to generate a 3D image of 4688 × 5160 × 4800 voxels. The testing was performed on a sub-volume of this image. Figures 4 and 5 show ESRGAN results of high perceptual quality that captured both rock textures and clays that were unresolved in the low-resolution image. While ESR-GAN results look almost virtually indistinguishable from their HR counterparts, we note that this study does not assess or report standard image-based metrics like PSNR and SSIM. As shown by Ahuja et al. (2022), such metrics can be misleading for micro-CT images trained using independent acquisitions and are not necessarily the most relevant metrics for the given context. For the same reasons, we too opted to evaluate our  For image segmentation, we used the Trainable WEKA Segmentation (Arganda-Carreras et al. 2017) that implements random forest (RF) machine learning algorithm. Berg et al. (2018) showed that with experienced training, trainable WEKA segmentation performed the best amongst seven different image-processing pipelines without any prior filtering. Benchmark study by Reinhardt et al. (2022) further showed that, compared to other machine learning and conventional segmentation techniques, user-bias of RF approaches was minimized due to constant interaction between experienced users and the RF classifier. However, the high classification uncertainty at phase transition voxels and tight pore connections due to partial volume effect can lead to increased user bias (Reinhardt et al. 2022). To take this uncertainty into account, we performed three realistic segmentations that an experienced user would be likely to consider. The segmentation "base" case in our study refers to the "ideal", high-confidence pore segmentation that would have been chosen if only a single segmentation was performed by the user for each image. With the base case as a starting point, two other cases "min" (porosity lower than "base") and "max" (porosity higher than "base") were trained by focusing on the pore-grain transition voxels (as shown in Figs. 6 and 7). This allowed us to assess a range of feasible segmentation results instead of relying on a single solution and further highlighted the higher susceptibility of LR images to user-bias for the same given volumes. Once the images are segmented, we compared the porosity and absolute permeability (computations performed with OpenFOAM) of the test subset. We observe that enhanced resolution of ESRGAN images reduced segmentation uncertainty (lower dispersion in porosity and permeability results) as we can see in Figs. 8 and 9. Although the base segmentations of the low-resolution images did not give bad results all the time, we notice that a different user can have segmented the images in a different way especially that the uncertainty range is very large for low-resolution images porosity and permeability estimations. We also notice that there were no major changes in the percolation thresholds between lowand super-resolution simulations. However, as the low-resolution percolation thresholds  were not very different from the high-resolution ones, we did not observe major changes on this parameter from the super-resolution images. This is not surprising as oil would invade large pores first during primary drainage and these are reasonably well resolved on all the images. Next, we check the bodies and throat size distributions for high-resolution, low-resolution, and super-resolution images. For that, we perform a pore network extraction using the maximum ball algorithm (Raeini et al. 2017) in order to determine the radii of the bodies and throats. We clearly observe that the estimation of pore and throat radii is considerably improved with super-resolution (Figs. 10 and 11).
Having validated the method on single-phase properties, we verified whether ESRGANgenerated images produce similar results to high-resolution images of the same volume on multiphase flow results. Therefore, we performed multiphase flow simulations using TotalEnergies' DRP simulation workflow that is based on pore network modelling technique (Regaieg et al. 2022(Regaieg et al. , 2021Regaieg and Moncorgé 2017). Figure 12 presents simulated primary drainage capillary pressure curves using low-resolution, high-resolution and super-resolution images for both Bentheimer and Reservoir C. As mentioned previously, we have tried to use several realistic segmentations for each case to illustrate the Fig. 9 Porosity (a), permeability (b) and percolation threshold (c) evolution depending on image resolution and segmentation hypothesis on Reservoir C images for factor 4 enhancement Fig. 10 Comparison between pore size distribution for bodies (up) and throats (down) using 3 images: low resolution (red), super-resolution (blue) and high resolution (green), and the three segmentation hypotheses (base, min and max) of pore network extracted from Bentheimer uncertainty related to this step, and therefore, we present envelopes of capillary pressure curves. We can clearly notice that ESRGAN makes our simulated primary drainage capillary pressure curves more accurate. We continue our analysis by comparing waterflooding relative permeability curves for several wettabilities. For simplicity, we only show the results of the base case segmentation; however, the conclusions do not change for the other segmentation scenarios. Figures 13 and 14 show water-wet and oil-wet relative permeability curves for simulations performed on low-resolution, high-resolution and super-resolution images. We highlight that we are comparing different wettability scenarios to make sure that our simulations are more accurate for different scenarios with different physics involved. In fact, characterization of wettability input is one of the major limitations of the PNM technique and is beyond the scope of this paper. Our proposed solution to tackle this challenging Fig. 11 Comparison between pore size distribution for bodies (up) and throats (down) using 3 images: low resolution (red), super resolution (blue) and high resolution (green), and the three segmentation hypotheses (base, min and max) of pore network extracted from Reservoir C Fig. 12 Comparison between computed capillary pressure curve envelopes for all three segmentation scenarios of the three datasets: low resolution (red), super resolution (blue) and high resolution (green) for: Bentheimer rock (a) and Reservoir Rock C (b) topic is addressed in Regaieg et al. (2022). For the water-wet case, we clearly observe that the residual oil saturation estimation improved considerably with the super-resolved image compared to the low-resolution one. This could be explained by the fact that snap-off capillary pressure thresholds are directly linked to the sizes of the throats that are not well captured in the low-resolution image. This is expected as in water-wet waterflood, the residual Fig. 13 Comparison between computed waterflood oil-wet (up) and water-wet (down) relative permeability curves on Bentheimer pore network using 3 images with base segmentation: low resolution (red), super resolution with a factor 4 (blue) and high resolution (green) in linear (left) and semi-log (right) Fig. 14 Comparison between computed waterflood oil-wet (up) and water-wet (down) relative permeability curves on Reservoir Sandstone C pore network using 3 images with base segmentation: low resolution (red), super resolution with a factor 4 (blue) and high resolution (green) in linear (left) and semi-log (right) oil saturation is generally the result of the competition between cooperative pore filling and snap-off. Underestimating snap-off events would result on lower S or as observed in our simulations. In the oil-wet case, the differences were less noticeable between the three simulations and residual oil saturation estimations were close in all scenarios. In this case, the S or is mainly controlled by the collapse of oil layers which is directly linked to the position of the water layers at the end of primary drainage, therefore, less related to the geometry, and this could explain why the results do not change notably for this case.
These encouraging results confirm the potential of the method, and therefore, we propose to apply it and test it for resolution improvement with a factor 8 in each direction.

Factor 8 Enhancement
First, we trained the neutral network with high-and low-resolution images of the same volume and from two different acquisitions of a Bentheimer sample with, respectively, 1.45 microns and 11.6 microns resolutions. We used 1000 slices to have 3000 crops for the training and 1000 crops for the cross-validation. Then, we used different 1000 slices for testing the accuracy of the method. Again, the training is performed using crops of 384 × 384 pixels with 36 GPUs and lasts 2 days. The super-resolved images were almost undistinguishable from high-resolution images and texture was well captured as we can see in Fig. 15. This is another example of the high perceptual fidelity achieved by ESRGAN.
To test the method, we performed single-and multi-phase flow simulations. However, as mentioned earlier, image segmentation is an uncertain process, and to take this uncertainty into account, we performed three realistic segmentations that an experienced user would be likely to consider. The segmentation philosophy is the same as that described above for factor 4 cases. The different segmentation realizations are shown in Fig. 16. Because of the coarser voxel size of LR image in × 8 case, we illustrate the increasing challenge of pore-grain boundary identification with a 2D zoom on the 3D results (Fig. 16), and its significant effect on the phase classification is reflected in the high porosity discrepancy of LR images in Fig. 17. Indeed, partial volume and resolution effects at such boundaries are the most significant contributors to DRP image processing sensitivity.
Next, we compared the porosity and pore size distributions in low-resolution, high-resolution and super-resolution images. We can clearly see there is less dispersion in the porosity estimation in the ESRGAN images compared to the low-resolution images (Fig. 17).

Fig. 15
Comparison between low-resolution (a), high-resolution (b) and super-resolution (c) cropped images of Bentheimer rock with a factor 8 resolution enhancement on the test dataset The pore size estimation and the average throat and body radii estimations improved drastically with super-resolution enhancement (Fig. 18). Additionally, we also performed singlephase flow simulations to compare the absolute permeability computed on low-resolution, Fig. 16 Example of the several segmentation hypotheses considered in this work for low resolution (LR), high resolution (HR) and super-resolution (ESRGAN) with factor 8 resolution enhancement on Bentheimer image Fig. 17 Comparison of permeability for low-, high-and super-resolution images for several realistic segmentations (max and min segmentation for each case) high-resolution and super-resolution images. We observe that the dispersion is reduced with ESRGAN and that the computed values are closer to the ones obtained from the highresolution images (Fig. 17). These observations are consistent with the conclusions of the previous section.
Subsequently, we perform multiphase flow simulations. Figure 19 shows the simulated primary drainage capillary pressure curves with the three images. As this curve is directly linked to the size of the pore bodies and the throats and to the connectivity of the sample, it is a good test for our method. It shows that the computations on the low-resolution image underestimated the capillary pressure especially for low water saturations. This is consistent with pore size distribution histograms described above where in the low-resolution images the pore sizes were overestimated. This issue was solved with ESRGAN and the simulated capillary pressure curve was in a good agreement with the calculations done on the high-resolution image. As in the simulations with factor 4 resolution enhancement, waterflooding relative permeability computations became more accurate with ESRGAN. Fig. 18 Comparison between pore size distribution for bodies (up) and throats (down) on Bentheimer pore network using 3 images: low resolution (red), super-resolution with a factor 8 (blue) and high resolution (green), and the three segmentation hypotheses (base, min and max) Fig. 19 Comparison between computed capillary pressure curves using 3 images with several realistic segmentation: low resolution (red), super resolution (blue) and high resolution (green) We could clearly observe that residual oil saturation estimations improved considerably for the water-wet case and a small improvement was observed in the computations of relative permeability for the oil-wet case (Fig. 20).

Testing and Validation of Pore Network Stitching
Super-resolution generates large images that are challenging to handle by the pore network extraction. Therefore, a stitching strategy has been developed as mentioned earlier in order to enable the extraction of a pore network from these large images and to enable performing large-scale flow simulations. To validate the stitching approach, the same sensitivity study on both the overlapping and transition length was performed on images of four rock samples (three reservoir sandstones (Fig. 21) and a Bentheimer outcrop).
As illustrated in Fig. 22, the validation process consisted in extracting two networks on the same image. The first was extracted using the full image and the second using the two subnetworks stitched together with the approach described previously.
It appears in all four cases that the overlapping length only had little impact on the obtained results compared to the transition length, as confirmed with the low dispersion on absolute permeability observed in Fig. 23. Moreover, a plateau in discrepancy shows the optimal cost of the stitching method.
When the transition length is too small (lower than roughly five times the mean throat length), the absolute permeability of the stitched network is lower than expected. This can be explained by the lack of connectivity in the transition length where long throats, which bypass the transition zone, are not considered during the stitching. When the transition zone tends to be too large, a small overestimation (less than 4%) of the absolute permeability is observed. For larger transition lengths, more inconsistencies are corrected in the transition zone, which leads to larger error in absolute permeability.

Fig. 20
Comparison between computed waterflood oil-wet (up) and water-wet (down) relative permeability curves using 3 images: low resolution (red), super resolution (blue) and high resolution (green) in linear (left) and semi-log (right) This is a result of the inconsistency of the extraction code (Raeini et al 2017) that we are using by cropping. These inconsistencies are not due to the method but only to the implementation that uses direction-dependent algorithm to read the 3D images. Therefore, these inconsistencies will cumulate as the transition zone gets larger. We are reimplementing Raeini et al (2017) code to remove the inconsistency by cropping, and this problem should be solved in the future. However, this error remains acceptable for all cases given in Fig. 23. Therefore, we consider a transition length of five times the mean throat length in the pore network stitches presented in this paper.
A similar comparison was performed on relative permeability on the same four rocks in both oil-wet and water-wet conditions. Figures 24 and 25 show the corresponding  and Reservoir Sandstone C (d). The different grey curves refer to different overlap and transition length combinations relative permeability curves; the different grey curves refer to overlap and transition length combinations considered in the sensitivity study. It allows to evaluate the maximum impact of stitching method on relative permeability compared to the reference. In each situation, the relative permeability in the stitched networks agreed well with the reference network.
Having tested a single stitching configuration does not guarantee that the addition of the error for several stitched networks remains reasonable. Therefore, a study of the consistency of the stitching relative permeability was conducted to ensure that the error remains reasonable even when the error of several stitching accumulates. A single extraction was performed on an image of Reservoir C rock with 935 × 935 × 4678 voxels; then, several stitching decompositions were performed and relative permeability curves were compared against the reference. Figure 26 shows a very good agreement between primary drainage relative permeability curves for the different stitched networks and the reference curves. As primary drainage follows an invasion percolation pattern, it is a good indication that the main connections are well conserved during the stitching and that cumulated error remains reasonable. Similarly, waterflood relative permeability curves are compared for water-wet and oil-wet scenarios (Fig. 27), and  Comparison between computed waterflood oil-wet (up) and water-wet (down) curves (in linear (a) and semi-log scale (b)) for several stitching decompositions (in 3D) with respect to reference single extraction Fig. 28 Comparison between computed primary drainage capillary pressure curves for several stitching decompositions (in 3D) with respect to reference single extraction overall, very good agreement was obtained between the stitched networks and our reference network. Figure 28 shows the primary drainage capillary pressure curves obtained from the stitched networks were very close comparing to the reference case. This gives us confidence about this approach.

Generation of Large Super-Resolved Images and Extraction of PNM
Having validated the ESRGAN and the stitching methodology, we applied them to generate large images with good resolution and extract a pore network from them. First, a low-resolution Bentheimer image at 5.8 microns was augmented using ESR-GAN to improve its resolution by a factor 4 in each direction and obtain an image with 3184 × 3280 × 12168 voxels (roughly 4800 3 voxels). Then, a stitching process was applied on the super-resolution image described previously after slightly cropping of the edges. The image was divided into several sub-volumes ( 2 × 2 × 9 ) along X, Y and Z, and networks were extracted from each sub-image and stitched together. A network with 3.6 million elements was generated (Fig. 29). Similarly, a low-resolution Reservoir Rock C image at 5.36 microns resolution was enhanced with ESRGAN to increase its resolution to 1.34 microns and generate an image with 4688 × 5160 × 4800 voxels (roughly 4900 3 voxels). Then, a stitching process was applied after dividing the image into sub-volumes of ( 4 × 3 × 4 ) along X, Y and Z, and networks were extracted from each sub-image and then stitched together to form a large pore network with 9.2 million elements (Fig. 29).
These newly generated pore networks are used in primary drainage pore network simulations (Fig. 30) where we present relative permeability curves for primary drainage on Bentheimer and Reservoir Rock C. We can notice that the flow for the sample with higher aspect ratio had higher percolation threshold comparing to the simulations on the cubic volume. We also plot the primary drainage capillary pressure curves for both rocks (Fig. 31) and observe that it was higher for the Reservoir Rock C which is expected as it had smaller pores than Bentheimer.

Conclusions
In this work, we applied ESRGAN to micro-CT images of rocks with resolution improvement factors of ×4 and ×8 in each direction. The method was tested on two rocks and it provided visually indistinguishable results with very high perceptual quality and ESRGAN was able to reproduce image rock textures and small pore features of previously unresolved areas (e.g. clays) that were similar to the original HR images. This helped us to reduce the segmentation uncertainty, and we observed lower dispersion in porosity and permeability results than in low-resolution images. We also observed that ESRGAN improved the pore geometry characterization and drainage Pc and pore-size-distribution results showed very good agreement with HR results compared to LR results.
Multiphase flow simulations were performed on high-resolution and low-resolution images and showed that capillary pressure and relative permeability computations were more accurate with ESRGAN. Even for a relatively homogeneous and well-connected rock like Bentheimer and Reservoir Sandstone C, the improvements in end-points and residual oil saturation value for relative permeability results of super-resolution image (compared to LR results) are significant. Using ESRGAN, large and highly resolved 3D images were generated for Bentheimer and Reservoir Rock C with sizes of 3184 × 3280 × 12168 and 4688 × 5160 × 4800 voxels, respectively. In order to enable the extraction of a pore network from a large image, a new stitching process was developed on networks extracted from overlapped sub-volume of a given image. We tested this method on single-and multiphase flow simulations for a single or several stiches, and we observed that even with large number of stitchings the computed properties were still accurate.
Finally, we used the stitching methodology to extract large pore networks from the super-resolved images and we obtained pore networks with, respectively, 3.6 and 9.2 million elements. These networks were used to perform primary drainage multiphase flow simulations. This work enables us to extract pore networks that we can use in large multiphase flow simulations. However, before that we need to solve the other weak point of PNM, which is wettability characterization (Regaieg et al. 2022).