Zhongzheng Wang
email: zhongzheng.wang@qut.edu.au

Hyogu Jeong

Yixiang Gan

Jean-Michel Pereira

Yuantong Gu
email: yuantong.gu@qut.edu.auc

Emilie Sauret
email: emilie.sauret@qut.edu.au

Pore-scale modeling of multiphase flow in porous media using a conditional generative adversarial network (cGAN)

Multiphase flow in porous media is involved in various natural and industrial applications including water infiltration into soils, carbon geosequestration, and underground hydrogen storage. Understanding the invasion morphology at pore scale is critical for better prediction of flow properties at the continuum scale in partially saturated permeable media. The deep learning method, as a promising technique to estimate the flow transport processes in porous media, has gained significant attention. However, existing works have mainly focused on single-phase flow, whereas the capability of data-driven techniques has yet to be applied to the pore-scale modeling of fluid-fluid displacement in porous media. Here, the conditional generative adversarial network (cGAN) is applied for pore-scale modeling of multiphase flow in two-dimensional porous media. The network is trained based on a data set of porous media generated using a particle-deposition method, with the corresponding invasion morphologies after the displacement processes calculated using a recently developed interface tracking algorithm. The results demonstrate the capability of data-driven techniques in predicting both fluid saturation and spatial distribution.

It is also shown that the method can be generalized to estimate fluid distribution under different wetting conditions and particle shapes. This work represents the first effort on the application of the deep learning method for pore-scale modeling of immiscible fluid displacement and highlights the strength of data-driven techniques for surrogate modeling of multiphase flow in porous media.

I. INTRODUCTION

Multiphase flow in porous media has been investigated extensively mainly in the context of geological systems for applications including enhanced oil recovery, carbon geosequestration, and underground hydrogen storage [1][2][3][4] . It has been established that the fluid-fluid displacement processes depend on both flow conditions and fluid properties, which was revealed in the phase diagram proposed in the seminal work by Lenormand 5 , where three regimes including capillary fingering, viscous fingering, and stable displacement were found to be governed by the viscosity ratio of the two fluids and Capillary number that reflects the relative importance of viscous force to capillary force. Besides, the significant impacts from the wettability of the porous media, quantified by the contact angles measured at threephase contact lines, have been identified experimentally and numerically [6][7][8][9][START_REF] Primkulov | Wettability and lenormand's diagram[END_REF] . Generally, as the wetting conditions of porous media shift from non-wetting (drainage processes with contact angle measured within the invading phase being less than 90 • ) to wetting (imbibition processes), smoother invasion fronts with less trapping of defending phase are observed, which is explained by the favored cooperative pore-filling events at pore scale [START_REF] Cieplak | Dynamical transition in quasistatic fluid invasion in porous media[END_REF][START_REF] Holtzman | Wettability stabilizes fluid invasion into porous media via nonlocal, cooperative pore filling[END_REF] .

The knowledge of multiphase flow in porous geological systems, i.e., rocks and soils, has been employed to study other topics involving different types of porous materials, e.g., gas diffusion layer in fuel cells [START_REF] Chapuis | Twophase flow and evaporation in model fibrous media: Application to the gas diffusion layer of pem fuel cells[END_REF][START_REF] Bao | Roughness effects of gas diffusion layers on droplet dynamics in pemfc flow channels[END_REF] . In most applications mentioned above, inferring macroscopic metrics of fluid-fluid displacement processes, such as displacement efficiency, fractal dimension of the invading fluid distribution, and relative permeability from fluid properties and flow conditions has been a central focus. In the past decades, increased emphasis has been put on the deterministic prediction and control of fluid transport phenomena especially in artificial porous media where the geometry of pore structure can be controlled, such as microfluidic devices for applications including drug delivery and selective metallization [START_REF] Zhao | Multiphase flow microfluidics for the production of single or multiple emulsions for drug delivery[END_REF][START_REF] Dudukovic | Cellular fluidics[END_REF] . Specifically, how geometry and topology of the pore structure affects multiphase flow remains an active area of research [START_REF] Holtzman | Effects of pore-scale disorder on fluid displacement in partially-wettable porous media[END_REF][START_REF] Rabbani | Pore geometry control of apparent wetting in porous media[END_REF][START_REF] Wang | Disorder characterization of porous media and its effect on fluid displacement[END_REF][START_REF] Hu | Phase diagram of quasi-static immiscible displacement in disordered porous media[END_REF][START_REF] Shi | Effects of topological disorder in unsaturated granular media via a pore-scale lattice Boltzmann investigation[END_REF][START_REF] Wang | Effect of grain shape on quasi-static fluid-fluid displacement in porous media[END_REF][23] . Recent works have shown that, by careful design of porous media with a particle size gradient, the capillary or viscous fingering can be suppressed during fluid-fluid displacement processes [START_REF] Rabbani | Suppressing viscous fingering in structured porous media[END_REF][START_REF] Lu | Controlling capillary fingering using pore size gradients in disordered media[END_REF] . Despite these efforts, the deterministic prediction of multiphase flow in porous media remains elusive.

Various numerical approaches have been employed to supplement experiments in understanding the fluid transport in porous media at pore-scale, including conventional Navier-Stokes equations solvers such as volume of fluid method [START_REF] Huang | Computer simulation of two-phase immiscible fluid motion in unsaturated complex fractures using a volume of fluid method[END_REF][START_REF] Raeini | Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method[END_REF][START_REF] Ambekar | Pore-resolved volume-of-fluid simulations of two-phase flow in porous media: Pore-scale flow mechanisms and regime map[END_REF] , mesoscale methods such as the lattice Boltzmann method [START_REF] Kruger | LBM -the principles and methods[END_REF][START_REF] Chen | A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications[END_REF][START_REF] Perumal | A review on the development of lattice Boltzmann computation of macro fluid flows and heat transfer[END_REF] , and pore-network models [START_REF] Cieplak | Dynamical transition in quasistatic fluid invasion in porous media[END_REF][START_REF] Blunt | Relative permeabilities from two-and three-dimensional pore-scale network modelling[END_REF][START_REF] Valvatne | Predictive pore-scale modeling of two-phase flow in mixed wet media[END_REF] . On the other hand, with the recent drastic advancement in the development of data-driven techniques, deep learning has become a promising tool for efficiently modeling fluid transport in porous media, where focus was placed on estimating the permeability of single-phase flow through a porous medium using artificial neural networks (ANNs) [START_REF] Rabbani | Hybrid pore-network and lattice-boltzmann permeability modelling accelerated by machine learning[END_REF] , convolutional neural networks (CNNs) [35][START_REF] Sudakov | Driving digital rock towards machine learning: Predicting permeability with gradient boosting and deep neural networks[END_REF][START_REF] Alqahtani | Machine learning for predicting properties of porous media from 2d x-ray images[END_REF][START_REF] Hong | Rapid estimation of permeability from digital rock using 3d convolutional neural network[END_REF][START_REF] Rabbani | Deepore: A deep learning workflow for rapid and comprehensive characterization of porous materials[END_REF][START_REF] Alqahtani | Flow-Based Characterization of Digital Rock Images Using Deep Learning[END_REF][START_REF] Elmorsy | Generalizable permeability prediction of digital porous media via a novel multi-scale 3d convolutional neural network[END_REF][START_REF] Wang | A physics-informed and hierarchically regularized data-driven model for predicting fluid flow through porous media[END_REF][START_REF] Kamrava | Physics-and image-based prediction of fluid flow and transport in complex porous membranes and materials by deep learning[END_REF][START_REF] Jin | Prediction model of velocity field around circular cylinder over various reynolds numbers by fusion convolutional neural networks based on pressure on the cylinder[END_REF][START_REF] Pawar | Data-driven recovery of hidden physics in reduced order modeling of fluid flows[END_REF] , or hybrid methods [START_REF] Kamrava | Linking morphology of porous media to their macroscopic permeability by deep learning[END_REF][START_REF] Tian | Permeability prediction of porous media using a combination of computational fluid dynamics and hybrid machine learning methods[END_REF] . To further obtain the spatial flow details, approaches based on CNNs have been proposed to either accelerate or directly predict the velocity field [START_REF] Hennigh | Lat-net: Compressing lattice boltzmann flow simulations using deep neural networks[END_REF][START_REF] Wang | Ml-lbm: Machine learning aided flow simulation in porous media[END_REF][START_REF] Ribeiro | Deepcfd: Efficient steadystate laminar flow approximation with deep convolutional neural networks[END_REF] . A recent review on the development of deep learning techniques on pore-scale modeling is given by Wang et al. [START_REF] Wang | Deep learning in pore scale imaging and modeling[END_REF] . Despite tremendous progress in the applications of data-driven methods in modeling transport problems in porous media, most studies have focused on single-phase flow. For multiphase flow, Ganti et al. [START_REF] Ganti | Data-driven surrogate modeling of multiphase flows using machine learning techniques[END_REF] presented the first attempt on applying data-driven methods for surrogate modeling of diesel jet injected into a quiescent nitrogen environment. Wen et al. [START_REF] Wen | Towards a predictor for CO2 plume migration using deep neural networks[END_REF] adopted CNNs as a computationally efficient substitute for predicting multiphase flow in the context of carbon dioxide storage. Shokouhi et al. [START_REF] Shokouhi | Physicsinformed deep learning for prediction of co2 storage site response[END_REF] presented a physics-informed deep learning method for prediction of CO 2 plume migration by modifying the loss function with the governing equations (continuity and Darcy's law). However, these aforementioned works are based on the equations at the continuum scale, neglecting the invasion mechanisms at pore scale. Nevertheless, it has been shown extensively that the pore-scale modeling of multiphase flow which includes the effects from wettability and pore structures is of vital importance for accurate prediction of fluid flow processes [START_REF] Rabbani | Pore geometry control of apparent wetting in porous media[END_REF][START_REF] Pavuluri | Towards pore network modelling of spontaneous imbibition: contact angle dependent invasion patterns and the occurrence of dynamic capillary barriers[END_REF][START_REF] Holtzman | Crossover from fingering to fracturing in deformable disordered media[END_REF][START_REF] Ju | Effects of pore characteristics on water-oil two-phase displacement in non-homogeneous pore structures: A pore-scale lattice Boltzmann model considering various fluid density ratios[END_REF][START_REF] Xu | Effect of pore geometry and interfacial tension on water-oil displacement efficiency in oil-wet microfluidic porous media analogs[END_REF][START_REF] Liu | Pore-scale modeling of spontaneous imbibition in porous media using the lattice boltzmann method[END_REF][START_REF] Reyssat | Imbibition in geometries with axial variations[END_REF][START_REF] Zacharoudiou | Pore-filling events in single junction micro-models with corresponding lattice Boltzmann simulations[END_REF][START_REF] Kou | Method for upscaling of co2 migration in 3d heterogeneous geological models[END_REF] .

In this work, the deep learning technique based on the conditional generative adversarial network (cGAN) called pix2pix by Isola et al. [START_REF] Isola | Image-to-image translation with conditional adversarial networks[END_REF] is applied for the pore-scale modeling multiphase flow in two-dimensional porous media. This method is chosen due to its wide applicability in handling image-to-image translation tasks. This is especially suitable for multiphase flow problems as both the input (porous structures) and the output (liquid distribution in the pore space) are often images such as those obtained using X-ray computed tomography. The flow condition considered in the work is the capillary-dominated regime, i.e., at low capillary number where the viscous effect is negligible compared with the interfacial tension, a scenario that can appear in applications such as carbon geosequestration and subsurface hydrogen storage [START_REF] Li | Morphology quantification of three-dimensional fluid invasion patterns[END_REF][START_REF] Hu | Experimental study on the displacement patterns and the phase diagram of immiscible fluid displacement in three-dimensional porous media[END_REF] . We firstly present the workflow, including generation of porous media, numerical method for fluid-fluid displacement, and the structure of cGAN. Next, the neural network is trained based on 4000 simulation results of drainage processes (porous media being non-wetting to the invading phase), and the prediction accuracy on test samples (containing 1000 simulation results) considering both statistical and deterministic metrics, i.e., saturation and invasion morphology, are examined. We explore the sensitivity of saturation and morphology prediction performance on the network architecture by varying the network depth. Then, we probe the generality of cGAN for multiphase flow by extending to cases with different wetting conditions and particle shapes. Finally, we discuss the implications of the findings and provide perspectives for potential future works.

overall framework

Porous Media Construction

II. METHODS

The workflow for predicting multiphase flow in porous media using deep learning techniques is shown in Fig. 1. Firstly, to obtain the required data set for training, the twodimensional porous media are generated as binary images, which are used as geometry input along with boundary conditions for multiphase flow simulation using numerical methods. The simulation results containing three phases, i.e., invading fluid, defending fluid, and solids are then binarized to represent the final morphology of the invading phase, which is regarded as the ground truth results. Note that, despite that both the defending phase and solids are in black color, the binary images of geometry input and ground truth result together should contain the complete information of phase distribution without loss of information, as both will be input into the learning machine, i.e., the conditional generative adversarial network (cGAN), during the training process. Once trained, given only a geometry input (test data) that the machine has not seen, the cGAN can make a prediction of invading phase distribution within the porous medium, which can be compared with the simulation results to examine the network performance. In the following, details are provided on the adopted methods in the workflow.

A. Porous Media Generation

Randomly arranged non-overlapping circular particles of different sizes have been commonly used as porous media with heterogeneous pore structures 9,66 . Here, we adopt a simple particle-deposition method to generate geometries with controlled statistical parameters.

Briefly, a circle with a prescribed radius R 0 is repetitively and randomly (uniform distribution) deposited into a square domain with an edge length L. The total number N of particles to be deposited can be calculated by N = (1-φ)L 2 πR 2 0 , with φ the porosity. During deposition, if the newly placed circle overlaps with an existing one, coalescence takes place, i.e., the particles merge into one larger particle with conserved total area. The location of the merged particle is calculated based on the area-weighted position, i.e.,

x new = (x 1 A 1 +x 2 A 2 )/(A 1 +A 2 ),
with x i the center coordinates and A i the area. A schematic showing the coalescence of two particles is shown in Fig. 2(A). A sample generated porous medium with a porosity φ = 0.6 is shown in Fig. 2(B). Interestingly, it is found that the resulting particle size distribution using this simple-deposition method can be well described by an exponential distribution [Fig. 2(C)]. To further control the geometrical features of porous media, firstly, a distance of inhibition h can be introduced to prevent particles from being too close, i.e., the coalescence criterion (ensuring non-overlapping) D < R 1 +R 2 is changed to D < R 1 +R 2 +h, with D the center-to-center distance between two particles, and R 1 and R 2 the particle radii. Note that h = 0 corresponds to non-overlapping criterion. The specification of an inhibition distance ensures sufficient pore spacing among particles such that the non-overlapping geometries can be accurately captured once the generated porous medium is binarized with a certain resolution. The other modification to control the particle size distribution is the introduction of a number of repetition N r during deposition to control the particle size distribution. Specifically, when depositing a new particle, if the coalescence criterion is triggered, coalescence is suppressed and another random location is chosen for the particle. This process repeats until a non-overlapping location is chosen. The number of repetition N r specifies the maximum times the "seek-new-location" process can be executed. Thus, N r = 1 corresponds to the original algorithm, and greater N r encourages space filling of small particles and avoids coalescence, which would lead to a more uniform particle size distribution. Fig. 2(D) shows three generated porous media in square domains of unit size with porosity φ = 0.6, R 0 = 0.03, h = 0.01, and N r = {1, 10, 100}, respectively. In this study, N r = 10 is chosen.

B. Interface Tracking Algorithm

To simulate multiphase flow in porous media, a recently developed interface tracking algorithm is used [START_REF] Wang | A pore-resolved interface tracking algorithm for simulating multiphase flow in arbitrarily structured porous media[END_REF] . The method is developed based on the algorithm originally proposed by Cieplak and Robbins 11 which considers essential pore-scale instability events, including burst, touch, and overlap. The original method has been successfully applied to reproduce experimental results of multiphase flow in porous media filled with perfectly spherical particles [START_REF] Hu | Phase diagram of quasi-static immiscible displacement in disordered porous media[END_REF][START_REF] Zhao | Capillary pinning and blunting of immiscible gravity currents in porous media[END_REF] . The key extension in the new algorithm is the consideration of unpin event, which results from the sharp edge pinning effect where the effective contact angle can be greater than the intrinsic one, and the contact angle at which unpin takes place can be determined according to a purely geometrical extension of Young-Dupre equation θ unpin = θ 0 + (180 • -α) with θ 0 the intrinsic contact angle and α the local corner angle [START_REF] Gibbs | The scientific papers[END_REF][START_REF] Oliver | Resistance to spreading of liquids by sharp edges[END_REF] . This is a common phenomenon observed in both natural and artificial surfaces where menisci can get pinned at sharp edges, e.g., Wu, Kharaghani, and Tsotsas 71 , Chen et al. [START_REF] Chen | Continuous directional water transport on the peristome surface of nepenthes alata[END_REF] . Fig. 3(A) shows the schematics of these pore-scale mechanisms, with light and dark blue curves representing the menisci position before and after an instability event, respectively. The main advantage of the interface tracking algorithm is its applicability to arbitrarily-structured porous media, as opposed to perfectly circular grains in the original model. This new method has been validated and applied to study fluid-fluid displacement processes in porous media with complex pore structures and different wetting conditions [START_REF] Wang | Effect of grain shape on quasi-static fluid-fluid displacement in porous media[END_REF][START_REF] Wang | A pore-resolved interface tracking algorithm for simulating multiphase flow in arbitrarily structured porous media[END_REF] . As illustrative examples, Fig. 3(B) and (C) show the simulated invasion morphologies of diagonal injection simulations in a representation of a Berea sandstone [START_REF] Guo | An experimental investigation of flow regimes in imbibition and drainage using a microfluidic platform[END_REF] with contact angles θ = {45 • , 165 • }, respectively. The algorithm is currently able to simulate multiphase flow in the capillarydominated regime, i.e., quasi-static processes with vanishing capillary number where the viscous effects are negligible. This regime is often encountered during fluid flow in fine soils and sands underground or in microfluidic devices where the size of pore space is small (in the order of micrometers). The algorithm is much more computationally efficient compared with conventional CFD methods, which facilitates the generation of data sets for training neural networks. A more detailed description of the algorithm can be found in the previous studies [START_REF] Wang | Effect of grain shape on quasi-static fluid-fluid displacement in porous media[END_REF][START_REF] Wang | A pore-resolved interface tracking algorithm for simulating multiphase flow in arbitrarily structured porous media[END_REF] and is omitted here for brevity.

C. Conditional Generative Adversarial Network

We adopt a conditional generative adversarial network (cGAN) called pix2pix as described in Isola et al. [START_REF] Isola | Image-to-image translation with conditional adversarial networks[END_REF] . This method is chosen due to its wide applicability in handling image-toimage translation tasks and ease of adoption without the need for parameter tweaking. This is especially suitable for multiphase flow problems as both the input (porous structures) and the output (liquid distribution in the pore space) are often images such as those obtained using X-ray computed tomography, after image treatment and segmentation (see e.g. Blunt et al. [START_REF] Blunt | Pore-scale imaging and modelling[END_REF] , Bruchon et al. [START_REF] Bruchon | Full 3d investigation and characterisation of capillary collapse of a loose unsaturated sand using x-ray ct[END_REF] ). A schematic of the cGAN is shown in Fig. 4, where the generator G, for a given input porous medium x, predicts the invading fluid distribution at percolation y, i.e., when the invading fluid reaches the outlet, while its adversary, the discriminator D, tries to classify whether the output result is "real" or "fake". In other words, during the training process, G is trained to produce output images that can "fool" D, whereas D is trained to distinguish fake images from real ones. Specifically, the architectures of the generator and discriminator are based on "ResU-net" and "PatchGAN" described in detail in Zhang, Liu, and Wang [START_REF] Zhang | Road extraction by deep residual u-net[END_REF] and Nie et al. [START_REF] Nie | TopologyGAN: Topology Optimization Using Generative Adversarial Networks Based on Physical Fields Over the Initial Domain[END_REF] , respectively. The exact network structures of the generator and the discriminator are provided in the Appendix (Fig. A1 and Fig. A2, respectively). All neural networks used in this framework are developed using Tensorflow [START_REF] Abadi | Tensorflow: A system for large-scale machine learning[END_REF] .

The loss function of the cGAN is given by:

L cGAN (G, D) = E x,y∼P data(x,y) [log D(x, y)] + E x∼P data(x) [log (1 -D(x, G(x)))], (1)
where G tries to minimize this objective against its adversary D that tries to maximize it.

At the same time, the generator not only learns to fool the discriminator but also is tasked to generate output close to the ground truth results in an L1 sense:

L L1 (G) = E x,y∼P data(x,y) [||y -G(x)|| 1 ]. ( 2
)
L1 norm rather than L2 is used to encourage clear boundary for the generator output [START_REF] Isola | Image-to-image translation with conditional adversarial networks[END_REF] , i.e., a clear boundary between fluid-fluid and fluid-solid interfaces. Hence, the final objective is [START_REF] Isola | Image-to-image translation with conditional adversarial networks[END_REF] :

G * = arg min G max D L cGAN (G, D) + λL L1 (G). (3)
In the original work of Isola et al. [START_REF] Isola | Image-to-image translation with conditional adversarial networks[END_REF] , it is shown that the quality of prediction is affected by the choice of λ. Particularly, L1 alone, corresponding to λ being too large or the absence of the first term on the right-hand side of Eqn.

(3), leads to blurry results, which is undesirable in the current application of immiscible fluids; A small value of λ, on the other hand,

gives sharper results but could introduce unwanted artifacts. Thus, a λ = 100 is chosen to reduce both of these artifacts, the same value adopted in Isola et al. [START_REF] Isola | Image-to-image translation with conditional adversarial networks[END_REF] . We believe that the optimal value of λ is problem-based (could be different depending on the application), which is worth further investigation. The standard approach is followed during training 79 : one gradient descent step is executed alternatively between D and G, and the objective is divided by 2 while optimizing D, which slows down the rate at which D learns compared with G. The Adam optimizer is adopted [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] , with a learning rate of 0.0002, and momentum parameters β 1 = 0.5, β 2 = 0.999. For the input, the generated porous media are binarized into 256x256 images as input images x for training. For flow simulation using the interface tracking algorithm, a diagonal injection setting is adopted, i.e., point inlet at the bottom left and point outlet at the top right, a typical geometry relevant to oil recovery processes.

The simulation results y are also binarized into 256x256 images where the morphologies of the invading fluid at percolation are colored in white. 5000 simulation cases under drainage conditions with a contact angle of 165 • are simulated, among which 4000 are for training, and 1000 for test.

III. RESULTS AND DISCUSSION

A. Performance Evaluation

A sample multiphase flow prediction for given input geometry from cGAN as well as the ground truth results from the simulation are displayed in Fig. 5(A). In the case of drainage processes (θ = 165 • ), the occupation of relatively big pores by the invading phase (white) is observed. Despite that there is some unphysical scattered invading phase that is disconnected from the main invading cluster, the prediction of phase distribution from cGANs compares well with the ground truth result. The image contrast in Fig. 5(A), i.e., cGAN prediction subtracted from the ground truth result, highlights the regions where the prediction and the ground truth result are different. Specifically, 1, 0, and -1 in the contrast map respectively represent the region of false invasion (not invaded in the ground truth), exact match, and false non-invasion (invaded in the ground truth). It can be seen that most of the unmatched regions (clusters in red or blue) are located at the menisci interfaces and are in the form of thin slices of the size of one or two layers of pixels. This, however, is not surprising as uncertainties are likely to be introduced when the geometry input and ground truth result are binarized into images with limited resolution. In the current work, with the choice of image resolution of 256x256 pixels, the characteristic geometrical feature size of the porous media, i.e., the average face-to-face distance between particles (or the average throat size), corresponds to 13.02 pixels. Higher image resolution associated with a greater number of pixels for the characteristic size may improve the accuracy at a compromise of increased computational cost.

The performance of cGAN in predicting multiphase flow is quantitatively evaluated based on two metrics. The first is the invading fluid saturation S inv = A inv /A pore , with A inv and A pore the area of invading fluid and total pore space, respectively. This is a critical parameter for the macroscopic characterization of multiphase flow in many engineering applications, such as enhanced oil recovery and carbon geosequestration. Note that, as the pore space is initially filled with only the defending phase, the value of S inv is also equal to the sweep efficiency of the displacement process, i.e., the proportion of the defending fluid that is displaced out of the porous domain. The comparison of saturation prediction of cGAN for the training data (blue) and test data (red) is shown in Fig. 5(B). The predictions on test data are clearly more sparsely distributed. However, a small standard deviation of 0.068 (black-dashed lines) of the error distribution shown in the inset indicates overall outstanding performance in saturation prediction. Indeed, the coefficient of determination is calculated to be R 2 = 0.981 on test data. The other metric, which is comparatively harsher compared with A inv , is the area ratio of total mismatched liquid distribution to the total pore space:

m = Ã A pore , (4)
with à the mismatched area [red and blue regions in the image contrast in Fig. 5(A)].

Therefore, m not only concerns the prediction of the total amount of invading liquid that is occupying the pore space after the displacement process, but also where it is distributed, i.e., the morphology, within the pore space. test data. The CDF shows the proportion of predictions (y-axis) that have errors that are less than a specific value (x-axis). A shift of curves towards the right indicates an increase in the error. It can be seen that the error in morphology prediction is higher than in saturation prediction. Specifically, if 20% is chosen as the error threshold, the corresponding percentage of saturation (P 20 s ) and morphology (P 20 m ) predictions are 98.7% and 89.4%, respectively. In other words, about 90% of predictions in the phase occupation status in the pore space have less than 20% error.

B. Sensitivity on Network Complexity

To investigate how the network architecture complexity impacts the prediction accuracy, the network depth of the cGAN generator is varied, which is illustrated in Fig. 6(A architecture-wise comparison, cGAN-16 is significantly outperformed by cGAN-4 and cGAN-8, whereas there is no significant difference between cGAN-4 and cGAN-8, indicating a saturated performance at cGAN-8, and further network complexity will not lead to a further increase in prediction accuracy. The quantitative metrics for performance quantification are summarised in Table . I. We note that apart from the architecture complexity examined in this section, other aspects of the neural network as well as parameters during training such as the type of activation function, the learning rate, and the momentum parameters could also impact the training time and the performance of the trained model; this is however beyond the scope of the current study.

C. Predicting Multiphase Flow in Porous Media

Since wettability and particle shape have been identified as key factors that can influence multiphase flow processes in porous media, it is important to assess the applicability and performance of cGAN in predicting multiphase flow under different wetting conditions and pore structures. So, we generate another 15000 cases of simulation results to cover different contact angles θ = {60 • , 105 • } and angular particles with sharp corners (represented by square shapes in this work). The porous media with angular grains can be generated based on the original circular ones by converting circular grains to square ones while preserving particle size (conserved area). Random rotation is then applied to every square particle.

Again, for each 5000 data sets, 4000 are for training and 1000 are for validation. cGAN for different wettability and particle shapes. Qualitatively, the predictions of the invading phase distribution within the pore space agree reasonably well with the simulation results. A general trend of more compact area occupied by the invading phase can be observed for cases with smaller contact angles, consistent with the expected transition from capillary fingering to stable displacement as contact angle decreases due to favored cooperative pore-filling events [START_REF] Cieplak | Dynamical transition in quasistatic fluid invasion in porous media[END_REF][START_REF] Holtzman | Wettability stabilizes fluid invasion into porous media via nonlocal, cooperative pore filling[END_REF] . Again, small regions of dispersed invading fluid that are disconnected from the main invading cluster can be observed. It is possible to design and carry out a post-processing procedure during which these isolated clusters are filtered out. However, we note that such a procedure may not be necessarily beneficial in interpreting the predicted results from cGAN. This is because the isolated clusters are not necessarily unphysical and might appear in three-dimensional displacement processes. In strong drainage, snap-off events caused by the swelling of defending phase corner flow could occur, which can disconnect the invading fluid [START_REF] Reynolds | Dynamic fluid connectivity during steady-state multiphase flow in a sandstone[END_REF][START_REF] Spurin | Mechanisms controlling fluid breakup and reconnection during two-phase flow in porous media[END_REF] , although the disconnection tends to be intermittent.

In strong imbibition, the fluid invasion process in the capillary-dominated regime follows a percolation-like pattern where the invading phase can be mainly transported through corner flows 9, [START_REF] Blunt | Multiphase Flow in Permeable Media: A Pore-Scale Perspective[END_REF] . This implies that the invading fluid could appear to be disconnected for given limited imaging resolution where the corner or thin film flows cannot be resolved. Therefore, in this work we refrain from applying any post-processing techniques to the results from predictions and the ground truth results. Specifically, the coefficients of determination are R 2 = {0.849, 0.956, 0.981} for θ = {60 • , 105 • , 165 • }, respectively, which also corresponds to P 20 s = {0.775, 0.940, 0.987} (proportion with less than 20% error in saturation prediction), respectively. One reason that may explain the decrease in accuracy as the porous media become more wetting to the invading fluid could be the favored non-local cooperative pore-filling events (or overlap event), since the pore invasion triggered by the overlap of two menisci depends on the exact location of both menisci, and as mentioned previously significant proportion of morphological inconsistency are located at the fluid-fluid interfaces [see the image contrast in Fig. 5(A)]. The average values of S inv for different contact angles are indicated as solid dots. As expected, the invading phase saturation increases as θ decreases, consistent with existing literature 8, [START_REF] Wang | Disorder characterization of porous media and its effect on fluid displacement[END_REF][START_REF] Geistlinger | The impact of wettability and surface roughness on fluid displacement and capillary trapping in 2d-and 3d-porous media part 1: Wettabilitycontrolled phase transition of trapping efficiency in glass beads packs[END_REF] .

The performance of cGAN on predicting fluid morphological distribution for different wettability and particle shapes is evaluated by plotting the CDF of error in invasion morphology predictions (Fig. 9). The overall trend of accuracy as a function of wettability is similar: greater error is observed when θ decreases. For circular grains with θ = 60 • , only 56% morphology predictions have less than 20% error, a noticeable reduction compared with the 91% for θ = 165 • . For porous media with angular particles, the proportion of predictions having less than 20% error is slightly smaller (86%) compared with circular grains at the same wetting condition (91%). This might be attributed to the decrease in the minimum throat size when converting circular grains into square ones. As mentioned previously, the porous media with square particles are generated by converting each circular grain to square ones with the area of each particle remaining the same (consequently the same porosity), after which a random rotation is applied to each individual particle. As a result, the minimum feature size of the porous media -the minimum possible face-to-face distance between two particles -is smaller for square particles compared with circular ones. Given the intrinsic uncertainties (inversely related to the image resolution) remain the same, slightly greater errors are observed in porous media with angular particles.

D. Discussion

The capability of data-driven method for predicting fluid-fluid displacement processes in porous media with different wetting conditions and particle shapes has been clearly demonstrated. The saturation estimation from cGAN is in good agreement with the ground truth results, whereas the prediction of spacial liquid distribution fails to achieve the same level of quantitative accuracy. The noticeable reduction in fluid morphology prediction highlights an important feature of multiphase flow in porous media compared with single-phase flow. During single-phase flow, small variations in local pore size is not expected to have a significant impact on the permeability. However, for multiphase flow, the global spatial liquid distribution can be very sensitive to local pore/throat size fluctuations, especially in processes where the capillary force and pore geometry governs the invasion sequence. An infinitesimal amount of variation in the size of one throat/pore could theoretically lead to the change of phase occupation status of an entire region, e.g., when all the subsequent invasion of a region depends on the filling of a specific pore/throat. An example of drastic change in the invasion morphology during imbibition processes due to slight variation in porous medium porosity or wettability has been reported recently [START_REF] Wang | Emergence of unstable invasion during imbibition in regular porous media[END_REF] . This, combined with the previously mentioned more demanding accuracy for the exact menisci location prediction during the cooperative pore-filling event, results in deteriorated performance in fluid morphology prediction during imbibition processes (θ = 60 • ). Despite the complexity involved in fluid-fluid displacement processes, the overall performance of cGAN on fluid morphology prediction is found to be satisfactory: 56% (imbibition) to 91% (drainage) correct prediction with 20% error tolerance, and 94% (imbibition) to 100% (drainage) with 40% error tolerance. In future works, images with enhanced quality, i.e., greater pixel density for the pore geometry, can be adopted to explore the improvement in pore-scale accuracy during imbibition processes.

In the current study, conditional generative adversarial networks were chosen as the network architecture for the pore-scale modeling of fluid-fluid displacement processes under different wettability and particle shapes. This is because, as mentioned before, (1) its proven strength of handling image-to-image translation tasks, and (2) the inputs and outputs of multiphase flow in porous media are often images as those obtained from X-ray CT.

By exploring the performance of cGAN under different network complexity, it is shown that the prediction accuracy is saturated at network cGAN-8, which is associated with a network depth of 6. However, we note that this observation is likely to be case-specific, i.e., the network depth needed for saturated performance depends on the smallest geometrical features of the pore structure, e.g., the characteristic throat size of the porous media. Here, cGAN has been focused on in this work, and we note that further studies are required for the exploration of optimal network structures using other alternative architectures. For the computational cost, the computation time is around 85 hours on an NVIDIA Tesla P100 GPU for the training process with 4000 training data. Once trained, the computation time for predicting 1000 test cases using cGAN is about 6 minutes, which is approximately two orders of magnitude faster than the interface tracking algorithm that was used for data generation [START_REF] Wang | A pore-resolved interface tracking algorithm for simulating multiphase flow in arbitrarily structured porous media[END_REF] . In the current work, despite that the same network architecture is adopted for all cases as shown in Fig. 7, the model is retrained for different wetting conditions and particle shapes. Future work could include embedding the contact angle into the training data such that the model can be used to predict the fluid distribution under various conditions. Alternatively, it is also possible to apply transfer learning to accelerate the training processes.

IV. CONCLUSIONS

We have presented a workflow where the data-driven technique is applied for predicting pore-scale multiphase flow in porous media, which is a crucial process in many engineering applications such as enhanced oil recovery and carbon geosequestration. It is shown that deep learning techniques can be used to not only predict the macroscopic metric during multiphase flow such as saturation, but also are able to perform satisfactorily in predicting spatial liquid distribution with different wetting conditions and particle shapes. Specifically, more than 99% (drainage) and 78% (imbibition) cases have less than 20% error for the saturation prediction, and more than 91% (drainage) and 56% (imbibition) cases have less than 20% error regarding the exact liquid spatial distribution. The greater error observed during imbibition processes is attributed to the more favored non-local cooperative porefilling events at smaller contact angles, a mechanism that is demanding in accurate prediction of the exact location of menisci. This work represents the first effort on the application of data-driven technique for the pore-scale modeling of immiscible fluid displacement in porous media. The direct consideration of wettability and pore structure geometry during multiphase flow distinguishes the current work from several recent studies in the literature [START_REF] Ganti | Data-driven surrogate modeling of multiphase flows using machine learning techniques[END_REF][START_REF] Wen | Towards a predictor for CO2 plume migration using deep neural networks[END_REF][START_REF] Shokouhi | Physicsinformed deep learning for prediction of co2 storage site response[END_REF][START_REF] Wen | U-fno-an enhanced fourier neural operator-based deep-learning model for multiphase flow[END_REF][START_REF] Yan | A robust deep learning workflow to predict multiphase flow behavior during geological co2 sequestration injection and post-injection periods[END_REF]

FIG. 1 .

 1 FIG. 1. Workflow of predicting multiphase flow in porous media using the conditional generative adversarial network (cGAN). Both the constructed porous media and simulation results (ground truth) of invasion morphologies are input into the cGAN for training. Once trained, given a new porous medium, cGAN directly predicts the invasion morphology of the invading phase at percolation.

FigureFIG. 2 .

 2 Figure 2 -Fractal circle generation

FigureFIG. 3 .

 3 Figure 3 -Interface Tracking Algorithm

FIG. 4 .

 4 FIG. 4. Conditional generative adversarial network (cGAN). The Generator G learns to generate the output that can fool the discriminator. The discriminator D learns to classify fake images (created by the generator) from the ground truth y.

FIG. 5 .

 5 FIG. 5. Error quantification for drainage with θ = 165 • . (A) From left to right: a sample input of a porous medium, the corresponding simulation result (ground truth), prediction from cGAN, and contrast between simulation result and cGAN prediction, respectively. (B) Comparison of invading phase saturation from simulation and cGAN prediction on training (blue) and test (red) samples. The coefficients of determination R 2 for training and test results are 0.998 and 0.981, respectively. The inset shows the distribution of relative error from cGAN on test samples, with the standard deviation 0.068 indicated by black-dashed lines. (C) Cumulative distribution function (CDF) for the error percentage for both saturation and morphology predictions. The percentage of cases that have less than 20% error on saturation-training, saturation-test, morphology-training, morphology-test are {100.0, 98.7, 99.9, 89.4} %, respectively.

  Fig. 5(C) shows the cumulative distribution function (CDF) of the errors in saturation ( s , absolute saturation difference between cGAN predictions and ground truth results) and morphology prediction ( m ) for both training and

FIG. 6 .

 6 Figure 6_v2 -network structure changes, and corresponding results

  ). The deepest bridging layers are progressively removed, which is associated with decreasing the network complexity of the generator. The obtained architectures are denoted as cGAN-4, cGAN-8, cGAN-16. Fig. 6(B) and (C) respectively show the error CFD for the saturation and morphology predictions with different network depths. As expected, compared with spatial fluid distribution, the saturation predictions are generally better for all cases. For

Fig. 7

 7 Fig.7shows three sample geometry inputs, ground truth results, and predictions from

  FIG. 8. Saturation prediction for multiphase flow under different wetting conditions θ = {60 • , 105 • , 165 • }, with R 2 = {0.849, 0.956, 0.981}, respectively. The dashed and solid curves surround the 80% and 50% cases based on the contours of density plot. The non-transparent circles denote the mean saturation for different θ.

  FIG. 9. Cumulative distribution function (CDF) for error percentage of invasion morphology predictions. For test cases, the proportion of cases that have less than 20% errors for circular grains θ = {60 • , 105 • , 165 • } and angular grains θ = 165 • are {0.56, 0.73, 0.91, 0.86}, respectively.

  FIG. A1. Schematic showing the network structure of the generator used in this work.

2

  M. Blunt, F. F.John, and F. M. Orr, "Carbon dioxide in enhanced oil recovery," Energy Conversion andManagement 34, 1197-1204 (1993), proceedings of the International Energy Agency Carbon Dioxide Disposal Symposium.

TABLE I .

 I Comparison of performance for the generator with different network depth. P 20 s and P 20 m represent the proportion of cases on the validation data set that have less than 20% error for saturation and morphology, respectively. s and m denote the average error for saturation and morphology predictions, respectively.

	Architecture	Network Depth	P 20 s	s	P 20 m	m
	cGAN-4	7	0.987	0.051	0.894	0.111
	cGAN-8	6	0.987	0.051	0.906	0.109
	cGAN-16	5	0.747	0.141	0.111	0.277

L. W. Lake, R. Johns, B. Rossen, and G. Pope, Fundamentals of Enhanced Oil Recovery (Society of Petroleum Engineers,

2014).

M. L. Szulczewski, C. W. MacMinn, H. J. Herzog, and R. Juanes, "Lifetime of carbon capture and storage as a climate-change mitigation technology," Proceedings of the National Academy of Sciences 109, 5185-5189 (2012).

Z. Kou, H. Wang, V. Alvarado, J. Fred McLaughlin, and S. A. Quillinan, "Impact of sub-core scale heterogeneity on co2/brine multiphase flow for geological carbon storage in

Real