Amir Hossein

Mohammadi Alamooti

Qumars Azizi

Hossein Davarzani
email: h.davarzani@brgm.fr.

A H Mohammadi Alamooti

Direct numerical simulation of trapped-phase recirculation at low capillary number

Keywords: Recirculation, volume of fluid, pore-doublet model, trapped phase, two-phase flow Highlights

come

Direct numerical simulation of trapped-phase recirculation at low capillary number

displacement mechanisms, fluid configurations, and dynamic properties. To simulate the porescale displacement at low capillary number, we used a filtered surface-force formulation of volume of fluid method, which was implemented in a separately available solver for OpenFoam package. The results showed that within the ranges of capillary number of invading phase analyzed in this study (in the order of 1×10 -7 to 1×10 -2 ), the recirculation phenomenon exists in trapped phases. During the imbibition mechanism, two stagnant regions are created adjacent to the fluid-fluid interface inside the invading fluid. Drag-force analysis on fluid-fluid interfaces shows that during imbibition the maximum force is exerted near the center of the interface, whereas during drainage more force is applied on two elongated interface tails on a solid surface. The centroids are elongated parallel to the interface during drainage and perpendicularly during imbibition, which is in concordance with drag-force distribution along

Introduction

Transport phenomena at pore-scale have gained great prominence in recent years since they are of significant importance in a variety of areas including enhanced oil recovery [START_REF] Zhao | Pore-scale modeling: Effects of wettability on waterflood oil recovery[END_REF], soil remediation [START_REF] Aminnaji | Effects of pore-scale heterogeneity on macroscopic NAPL dissolution efficiency: A two-scale numerical simulation study[END_REF], and carbon/hydrogen storage [START_REF] Ebigbo | A coupled, pore-scale model for methanogenic microbial activity in underground hydrogen storage[END_REF][START_REF] Riazi | Visualisation of mechanisms involved in CO2 injection and storage in hydrocarbon reservoirs and water-bearing aquifers[END_REF]. Various microscale forces including capillary, viscous, inertial, and gravitational forces play an important role in fluid dynamics especially at the interface [START_REF] Dejam | Reinfiltration through liquid bridges formed between two matrix blocks in fractured rocks[END_REF]Ferrari andLunati, 2014, 2013;[START_REF] Mashayekhizadeh | Pore-level observation of free gravity drainage of oil in fractured porous media[END_REF][START_REF] Méheust | Interface scaling in a twodimensional porous medium under combined viscous, gravity, and capillary effects[END_REF][START_REF] Raeini | Numerical modelling of sub-pore scale events in two-phase flow through porous media[END_REF]. The roles of rock and fluid properties that can significantly affect the displacement process are also important [START_REF] Avraam | Flow regimes and relative permeabilities during steadystate two-phase flow in porous media[END_REF]. Specifically, pore-scale investigation is a strong tool that can shed light on the subsurface process at the microscale. To investigate the physics of fluid flow at pore-scale, 2D glass micromodels can be used in conjunction with fluorescence microscopy and micro-particle image velocimetry (micro PIV) to directly observe and measure the pore-scale processes [START_REF] Jiménez-Martínez | Impact of saturation on dispersion and mixing in porous media: Photobleaching pulse injection experiments and shear-enhanced mixing model[END_REF][START_REF] Roman | Pore-scale visualization and characterization of viscous dissipation in porous media[END_REF][START_REF] Roman | Particle velocimetry analysis of immiscible two-phase flow in micromodels[END_REF]. [START_REF] Roman | Particle velocimetry analysis of immiscible two-phase flow in micromodels[END_REF] investigated the velocity profiles in two-phase flow through the heterogeneous sandstone-replica micromodel using micro PIV. The results depicted recirculation motion inside the immobile pockets of the wetting phase due to momentum transmitted by the flowing oil. [START_REF] Kazemifar | Quantifying the flow dynamics of supercritical CO2-water displacement in a 2D porous micromodel using fluorescent microscopy and microscopic PIV[END_REF] reported shear-induced circulation in trapped water ganglia during post-front passage displacement of water by supercritical CO2 in a homogenous micromodel. In concordance with the Kazemifar et al. work, Li et al (2017) observed the shear-induced circulation during post-front passage flow in a heterogeneous micromodel for a supercritical CO2/water system. They determined CO2 flow directions with the aim of circulation within the trapped water phase, as it is very difficult to directly measure CO2 flow. Using micro PIV, [START_REF] Heshmati | Interfacial boundary conditions and residual trapping: A porescale investigation of the effects of wetting phase flow rate and viscosity using microparticle image velocimetry[END_REF] observed the recirculation process in a pore-doublet model in the viscous dominant regime. They reported that the shear stress applied to the trapped non-wetting phase was intensified as the viscous forces increased. The effect of increasing injection velocity on the velocity of the non-wetting phase was less noticeable, indicating the presence of a slip boundary between the fluids. [START_REF] Zarikos | Velocity distributions in trapped and mobilized non-wetting phase ganglia in porous media[END_REF] analyzed the interaction of capillarity and momentum transfer between the two immiscible phases for different non-wetting trapped globules with various sizes. In contrast with stagnant areas where the viscous drag applied on the fluid-fluid interface is negligible, in nonzero velocity areas the momentum transfers cause circulations within the trapped phase.

The circulation velocity is maximum near the fluid-fluid interface and an increase in capillary number of invading phase increases the circulation velocity as well as the shift in circulation centroid toward high-velocity region. Similar to the circulation phenomenon process for trapped phases in porous media, there are other two-phase systems in which this process is noticeable [START_REF] Bhaga | Bubbles in viscous liquids: shapes, wakes and velocities[END_REF][START_REF] Maxworthy | Experiments on the rise of air bubbles in clean viscous liquids[END_REF]. The effect of capillary numbers, slug size, and viscosity ratios on internal recirculation within two phases in capillary channels have been investigated in several experimental [START_REF] Kashid | Liquid-liquid slug flow in a capillary: an alternative to suspended drop or film contactors[END_REF][START_REF] Liu | Micro-PIV investigation of the internal flow transitions inside droplets traveling in a rectangular microchannel[END_REF][START_REF] Ma | On the flow topology inside droplets moving in rectangular microchannels[END_REF] and numerical studies [START_REF] Kashid | Internal circulation within the liquid slugs of a liquid-liquid slug-flow capillary microreactor[END_REF][START_REF] Sarrazin | Experimental and numerical study of droplets hydrodynamics in microchannels[END_REF].

Contrary to the recent studies where the velocity fields were observed in only one phase (usually the water phase), two-phase flow can be explored using an appropriate two-phase flow simulator. Pore-network modeling as a low-cost computational tool has advanced during recent years in both quasi-static and dynamic modes [START_REF] Blunt | Flow in porous media-pore-network models and multiphase flow[END_REF][START_REF] Oren | Extending predictive capabilities to network models[END_REF]. Direct numerical simulation methods including smoothed particle hydrodynamic (SPH) as a meshless method [START_REF] Hirschler | Open boundary conditions for ISPH and their application to micro-flow[END_REF], the lattice Boltzmann method (LBM) [START_REF] Kang | Unified lattice Boltzmann method for flow in multiscale porous media[END_REF], and computational fluid dynamic methods (CFD) [START_REF] Raeini | Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method[END_REF] can perform simulation at high resolution of the real structure of pore spaces. [START_REF] Fakhari | A phase-field lattice Boltzmann model for simulating multiphase flows in porous media: Application and comparison to experiments of CO2 sequestration at pore scale[END_REF] using the lattice Boltzmann model for the simulation of multiphase flows of water/CO2 systems in a micromodel, observed the recirculating flows in the still-flowing CO2 phase as the impacts of inertial forces. Direct numerical simulation of pore-scale processes by CFD methods have gained great prominence due to their ability to simulate fluid flow for wide density and viscosity ratio ranges [START_REF] Meakin | Modeling and simulation of pore-scale multiphase fluid flow and reactive transport in fractured and porous media[END_REF]. In a recent study, using direct numerical simulation, the results of drainage processes in a 2D cavity showed better mixing processes and larger interfacial mass transfer in the presence of recirculation inside the cavity [START_REF] Maes | A new compressive scheme to simulate species transfer across fluid interfaces using the volume-of-fluid method[END_REF]. [START_REF] Shams | A study to investigate viscous coupling effects on the hydraulic conductance of fluid layers in two-phase flow at the pore level[END_REF] analyzed the viscous coupling effects as a function of the viscosity ratio, wettability, and varying fluid configurations in circular capillary tubes. [START_REF] Raeini | Numerical modelling of sub-pore scale events in two-phase flow through porous media[END_REF] studied the snap-off and layer flow to analyze the impact of geometry and flow rate on the hydraulic conductivity of disconnected ganglia using volume-of-fluid (VOF) based finite-volume method. To consider interfacial tension, various approaches can be applied to VOF; the continuous surface stress (CSS) method [START_REF] Gueyffier | Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF], the continuous surface force (CSF) method [START_REF] Brackbill | Dynamic immiscible displacement mechanisms in pore doublets: theory versus experiment[END_REF], and the sharp surface force (SSF) [START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF][START_REF] Raeini | Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method[END_REF]. In these approaches, it is very difficult to predict flow at low capillary numbers where the capillary forces are dominant, which can result in non-physical velocity currents around the interface [START_REF] Raeini | Numerical modelling of sub-pore scale events in two-phase flow through porous media[END_REF]. To eliminate the non-physical behavior from the numerical solution for the cases of capillary forces dominancy, the filtered surface force (FSF) is introduced [START_REF] Raeini | Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method[END_REF]. Specifically, the capillary pressure, dynamic pressure gradients, and viscous and inertial forces are separately calculated to prevent introducing instabilities in the numerical equations. It is shown that the spurious current is minimum around the interface when the FSF is used [START_REF] Pavuluri | Spontaneous imbibition in a microchannel: Analytical solution and assessment of volume of fluid formulations[END_REF].

The concept and application of the recirculation phenomenon are of great importance especially for processes contending with mass transport [START_REF] Kazemifar | Quantifying the flow dynamics of supercritical CO2-water displacement in a 2D porous micromodel using fluorescent microscopy and microscopic PIV[END_REF][START_REF] Li | Micro-PIV measurements of multiphase flow of water and liquid CO2 in 2-D heterogeneous porous micromodels[END_REF][START_REF] Maes | A new compressive scheme to simulate species transfer across fluid interfaces using the volume-of-fluid method[END_REF]. Recently some research has been carried out to investigate this phenomenon in both experimental [START_REF] Roman | Particle velocimetry analysis of immiscible two-phase flow in micromodels[END_REF][START_REF] Zarikos | Velocity distributions in trapped and mobilized non-wetting phase ganglia in porous media[END_REF] and numerical approaches [START_REF] Fakhari | A phase-field lattice Boltzmann model for simulating multiphase flows in porous media: Application and comparison to experiments of CO2 sequestration at pore scale[END_REF][START_REF] Maes | A new compressive scheme to simulate species transfer across fluid interfaces using the volume-of-fluid method[END_REF]. However, they have not discussed the role of displacement mechanism and fluid configurations on the recirculation process. Also, in the case of experimental investigations [START_REF] Kazemifar | Quantifying the flow dynamics of supercritical CO2-water displacement in a 2D porous micromodel using fluorescent microscopy and microscopic PIV[END_REF][START_REF] Roman | Pore-scale visualization and characterization of viscous dissipation in porous media[END_REF][START_REF] Roman | Particle velocimetry analysis of immiscible two-phase flow in micromodels[END_REF][START_REF] Zarikos | Velocity distributions in trapped and mobilized non-wetting phase ganglia in porous media[END_REF], the velocity profile is only presented inside the trapped phase. [START_REF] Maes | A new compressive scheme to simulate species transfer across fluid interfaces using the volume-of-fluid method[END_REF] investigated the role of recirculation phenomena on the mass transfer process for various viscosity ratios; however, they did not consider the effects of injection velocity and surface tension as other dynamic parameters involved in the capillary number. [START_REF] Zarikos | Velocity distributions in trapped and mobilized non-wetting phase ganglia in porous media[END_REF] investigated how the center of circulation and average velocity of the trapped phase moved for various capillary numbers, but they did not discuss the behavior of capillary number for various viscosity ratios. This work sheds light on the recirculation process for both imbibition and drainage mechanisms for various fluid configurations at low capillary number. Also, using drag force analysis on fluid-fluid and fluid-solid interfaces, we comprehensively discuss the behavior of recirculation phenomena. In addition, we discuss the analysis of dynamic parameters including injection velocity, viscosity ratio, and surface tension individually on the behavior of recirculation. We used the FSF formulation of the VOF method implemented in poreFoam (separately available solver for OpenFOAM) to simulate two-phase flow inside a 2D pore-doublet model at low capillary number. In this study, we addressed the question of how recirculation behavior changes under various dynamic and static conditions. As the multiphase-flow simulation is carried out in the real structure of pore spaces, this study is new to those that use simplified models such as pore network modeling. Also, using the FSF formulation of the VOF method minimizes the spurious currents introduced in numerical equations. On the other hand, the simulation run time is relatively long due to small time steps that should be used, especially at low capillary numbers. The main limitation is that even with this advanced method, the simulation of fluid flow for the invading phase capillary numbers lower than 10 -7 cannot be performed. Also, because of high computational cost, investigation of more complex porous media was not feasible. The increase in mesh size and shortness of time steps might have led to a long simulation time.

The results of this work can be applied for cases where recirculation plays a pivotal role along with mass transfer, such as CO2 storage. Analyzing drag force on the fluid-fluid interface is also important in processes where the balance of forces on the interface determines the mobilization or trapping of the droplet during enhanced oil recovery methods [START_REF] Dejam | Diffusive leakage of brine from aquifers during CO2 geological storage[END_REF][START_REF] Olayiwola | A comprehensive review on interaction of nanoparticles with low salinity water and surfactant for enhanced oil recovery in sandstone and carbonate reservoirs[END_REF][START_REF] Yang | Droplet mobilization at the walls of a microfluidic channel[END_REF][START_REF] Yang | Internal flow patterns of a droplet pinned to the hydrophobic surfaces of a confined microchannel using micro-PIV and VOF simulations[END_REF]. Figure 1

Numerical Method

Mathematical Formulation

This section described the mathematical model we used to solve two-phase flow. The Navier-Stokes equation describes the dynamics of two immiscible incompressible flow. Since part of the solutions are interface dynamics, a moving boundary problem is required [START_REF] Batchelor | An introduction to fluid dynamics[END_REF]. For an isothermal condition, the momentum-balance equation in terms of the viscous, inertial, capillary, and gravity forces is written as [START_REF] Raeini | Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method[END_REF]:

𝐷 𝐷𝑡 (𝜌𝑼) -𝜵. 𝑇 = -𝜵𝑃 𝑑 + 𝜌𝒈 + 𝒇 𝑐 -𝜵𝑃 𝑐 1
where U denotes the velocity field, 𝜌 the density, t the time, and g is gravity vector. The capillary pressure term is calculated using the equation:

𝜵. 𝜵𝑃 𝑐 = 𝜵. 𝒇 𝒄 3
With the boundary condition:

𝜕𝑃 𝑐 𝜕𝑛 𝑏 = 0 4
where 𝑛 𝑏 is the axis normal to the boundary. Equations ( 3) and (4) introduce a Neumann problem for capillary pressure term which has a unique solution. Including the effect of the capillary force in Navier-Stokes equations provides an opportunity to filter the instabilities introduced into the numerical equations where the capillary forces are high [START_REF] Raeini | Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method[END_REF].

At the fluid-fluid interface, the 𝛼 value is between zero and one. It is evolved with the advection equation of the form:

𝜕𝛼 𝜕𝑡 + 𝜵. (𝛼𝑼) = 0 5

By solving Equation 5, the indicator function is obtained then the capillary force 𝑓 𝑐 can be calculated as a body force [START_REF] Brackbill | Dynamic immiscible displacement mechanisms in pore doublets: theory versus experiment[END_REF]:

𝒇 𝑐 = 𝜎𝑘𝒏 𝑰 𝛿 𝑠 6
where 𝜎 is surface tension, nI is the unit interface normal vector, δ is the delta function active along the interface and 𝑘 = 𝛁. (𝒏 𝐼 ) is the interface curvature and is normal to the interface.

𝒏 𝐼 = 𝜵𝛼 |𝜵𝛼| 7
Detailed forced filtered formulation of the volume of fluid method is discussed in appendix A.

Computational Domain

Two-phase flow simulations were performed in designed 2D pore-doublets, inspired by the By subdividing the primary domain into non-overlapping subdomains as smaller discrete cells, flow equations 1-7 and A1-A11 are solved. The domain was meshed using OpenFOAM meshing tools snappyHexMesh [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF]. To visualize the processes, we used Paraview.

Fluid Properties

To analyze the role of displacement mechanisms on recirculation processes in various fluid configurations inside the pore assembly, we considered a set of fluid properties (Listed in Table 1). The properties were designed to ensure that enough shear stress was transmitted to the trapped phase in both drainage and imbibition mechanisms. To examine the effects of dynamic parameters involving capillary number including surface tension, velocity, and viscosity ratio, different ranges were considered (and explained in that section). As the Bond number (𝐵𝑜 = ∆𝜌𝑔𝐿 𝜎 ⁄ , where L is the characteristic length) is equal to zero due to the presence of fluids with the same densities, the effect of gravitational forces is not a focus in this study.

Boundary and Initial Condition

One of the essential physical parameters, which should be defined in the numerical solver, is suitable boundary conditions. Zero longitudinal velocity magnitude (Ux) is considered while the imposition of uniform velocity profile equal to 1.0×10 -4 m/s supplies the latitudinal velocity (Uy) at the inlet. For outlet boundary zero-gauge global pressure, we considered a resemblance with atmospheric pressure. Inlet gradient pressure values and the velocity gradient at the outlet were fixed at zero during all simulations. The grain surfaces predetermined wetting properties depending on the displacement mechanism i.e. fully-wet to the invading phase in the case of imbibition and fully non-wet to the invading phase in the case of drainage. The normal velocity component at a fluid-solid interface was set at zero (no-slip conditions). To analyze the recirculation behavior inside the trapped phases in similar conditions, the model was initially set to be filled with invading fluid then the trapped phase was locally patched at three individual locations on the right and left sides of the pore assembly and around the bottom dead-end. Then the invading fluid with a certain contact angle depending on the displacement mechanism (i.e.

zero for imbibition and 180º for drainage) was injected at the inlet.

Solution Methodology

For a time-dependent approach for 2D simulations, an open-source finite volume CFD package was used. The FSF formulation of the VOF method [START_REF] Raeini | Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method[END_REF] with pressure implicit-split operator (PISO) algorithms [START_REF] Issa | Solution of the implicitly discretised fluid flow equations by operator-splitting[END_REF] was a segregated procedure used to solve the pressure, velocity and volume fraction coupling equations. We used a semi-implicit formulation of capillary forces to confirm the accuracy and stability of the numerical method. Adjustable time-steps with a maximum courant number of 1.0×10 -5 [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF] were used.

The time steps were designed in a way to satisfy the constraints below [START_REF] Raeini | Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method[END_REF]:

𝛿𝑡 < 𝐶 𝛿𝑡 𝜇𝛿𝑥 𝜎 8
where 𝛿𝑡 is time-step size, 𝛿𝑥 is mesh resolution, and Cst is a constant varying depending on the formulation of capillary forces (explicit or semi-implicit) [START_REF] Raeini | Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method[END_REF]. The simulations in this work were performed on Intel Core i9, 5 GHz, 16 MB cache with 16 processors and 16 GB RAM.

Validation and Verification

Mesh Independency

The mesh density of the network is governed in a way that no dependency between the discretization and variable can be observed. The OpenFOAM snappyHexMesh utility was used for domain meshing. This tool allows good refinement of grids in the various parts of a pore assembly. The mesh quality was confirmed by analysis of grid refinement, where numbers of mesh elements were increased gradually from 8800 to 79,000, corresponding to the number of cells in the narrowest throat from 9 to 11, 13, and up to 23. The residuals, average velocity, and velocity distributions inside the sub-domains were calculated in each simulation to confirm the grid-size convergence. The results showed that the simulation with 19 cells in the narrowest throat equivalent to a mesh size of 55,635 cells led to minimum error while the finer mesh did not decrease the error. Figure 3 illustrates the mesh inside the pore assembly.

Results

This section considers the simulation results of two individual displacement mechanisms i.e.

imbibition and drainage through the pore-doublet model. In both cases three patches of defending phase were trapped in the right and left sides of the pore-doublet model as well as the bottom dead-end. This fluid configuration within the pore assembly not only facilitated the analysis of recirculation behavior in various arrangements to the flow direction but also provided the same conditions for both displacement mechanisms. Prior to any displacement simulation, the trapped patches were relaxed for 0.03 s to remove the capillary oscillation effects from the numerical calculations. During this process, no inflow or outflow of the fluids across the boundaries was considered, to make the interface reaching equilibrium configuration depending on the contact angle. Then the displacement process was simulated to reach the steady-state condition, depending on the displacement mechanism.

Model Validation

To ensure the accuracy of the numerical model, the simulation results were validated with available experimental data mentioned in the study of [START_REF] Heshmati | Interfacial boundary conditions and residual trapping: A porescale investigation of the effects of wetting phase flow rate and viscosity using microparticle image velocimetry[END_REF] [START_REF] Heshmati | Interfacial boundary conditions and residual trapping: A porescale investigation of the effects of wetting phase flow rate and viscosity using microparticle image velocimetry[END_REF]. They analyzed the recirculation process inside one side of a pore-doublet model (originally inspired by Chatzis and Dullien (1983)) for immiscible displacement. Table 2 shows the values for velocities at 10 µm intervals to the liquid-liquid interface inside invading and trapped phases. Our simulation results show that in a similar boundary and initial conditions to their experimental study, a recirculation zone was created inside the non-wetting trapped phase at the right side of pore assembly. Analysis of velocity distribution at both sides of liquid-liquid interface shows that the velocities are in good agreement with the experiment. Due to the low quality of images in [START_REF] Heshmati | Interfacial boundary conditions and residual trapping: A porescale investigation of the effects of wetting phase flow rate and viscosity using microparticle image velocimetry[END_REF], we validated this process using other experimental data. We used the experimental data related to a drainage process inside dead-end pore micromodels (i.e. eddy pockets) [START_REF] Roman | Pore-scale visualization and characterization of viscous dissipation in porous media[END_REF]. To simulate the drainage mechanism eddy pockets including three connected cavities were considered (Figure 4 (a)), then with the same boundary and initial conditions to their experimental study the non-wetting phase was injected through the micromodel. To avoid the borders effect, the central cavity results were selected to be compared with the experimental results. A mesh sensitivity was performed and mesh size equal to 93,600 cells (equivalent to cell size of 1.6 micron) was selected as the optimum density. As Figure 4 (b) shows, good agreement can be seen between the numerical simulation and experimental results, especially in phase distributions and the velocity profile inside the cavity. To better understand the velocity distribution inside the trapped phase, the values of Ux and Uy along the vertical axis to the centroid (dashed white line in Figure 4 (b,c))

for the simulation process are compared with the experimental result. The values of Uy are near zero for both experimental and simulation approaches. As can be seen in Figure 5 for Ux values, although the trends of experimental and simulation results are in concordance. There is only a minor change in the locations of the centroid (Ux=0) that can be attributed to the effect of differences in the degree of the sharpness of the corners in experimental and simulation analysis. Simulation Results (this study) elongated along the solid wall. Some research has been conducted to observe the characteristics of rotational flow during this mechanism in complex 2D porous media [START_REF] Fakhari | A phase-field lattice Boltzmann model for simulating multiphase flows in porous media: Application and comparison to experiments of CO2 sequestration at pore scale[END_REF][START_REF] Kazemifar | Quantifying the flow dynamics of supercritical CO2-water displacement in a 2D porous micromodel using fluorescent microscopy and microscopic PIV[END_REF][START_REF] Li | Micro-PIV measurements of multiphase flow of water and liquid CO2 in 2-D heterogeneous porous micromodels[END_REF]. Here the focus was on the various pore structures combined in a pore assembly to assess the recirculation process at low velocity. Using Equation 9, the capillary number of invading phase was calculated, where Ui is invading phase velocity, µi is invading phase viscosity, and σ is surface tension. The calculated 𝐶𝑎 𝑖 was around 3.3×10 - 5 and M, viscosity ratio, was 10.

Bottom Dead-End (behind the mean flow direction)

The presence of the bottom dead-end makes it possible to examine recirculation behavior inside a trapped phase in a location behind the mean flow direction. Figure 11 shows the fluid configuration and velocity fields in the bottom dead-end of the pore assembly. In contrast to two other cases, the dead-end is not fully trapped, and around half of it is occupied by the invading phase. This offers a chance to the invading phase to circulate inside the dead-end in both cases. Although in some studies more than one vortex is observed in a trapped phase [START_REF] Kashid | Internal circulation within the liquid slugs of a liquid-liquid slug-flow capillary microreactor[END_REF][START_REF] Liu | Micro-PIV investigation of the internal flow transitions inside droplets traveling in a rectangular microchannel[END_REF][START_REF] Zarikos | Velocity distributions in trapped and mobilized non-wetting phase ganglia in porous media[END_REF]

Drag-Force Analysis

Our goal is to quantify and analyze the velocity distributions inside the recirculating trapped phase in both imbibition and drainage mechanisms. Detailed velocity fields near the fluid-fluid interface are used for drag-force quantification [START_REF] Zarikos | Velocity distributions in trapped and mobilized non-wetting phase ganglia in porous media[END_REF]. The drag-force analysis helps us better understand recirculation behavior for various fluid configurations as well as the displacement mechanism. Viscous drag force per unit of the depth of interface exerted on each part of the interface is calculated as:

𝐷 = n. μ i ( ∂𝐮 𝐱 ∂y + ∂𝐮 𝐲 ∂x ) dl 10
Where D is drag force per unit of the depth of interface, n is the normal vector, μ i is the viscosity of the invading phase, and l is the interface length. The distribution of drag force per unit of interface depth for the right-side cavity is shown in Figure 12. It clearly shows that in the case of imbibition, the maximum drag force from the invading phase is exerted around the center of the interface, while in two stagnant regions at the upper and lower sides of the interface inside the invading phase, drag force has a reverse direction. It shows that in these areas the force is exerted from the recirculating trapped phase; however, it cannot cause recirculation inside the invading phase as the invading phase viscosity is noticeably high. In contrast, in the case of the drainage mechanism, drag force is at a minimum around the central parts of the interface and an increase can be seen by moving toward the upper and lower sides of the interface. This dragforce distribution along the interface directly affects the centroid shape inside the recirculating trapped phase. The values for the centroid elongation factor, which is defined as the parallel diameter of centroid to the interface divided by the perpendicular one, are 1.11 and 2.76 for imbibition and drainage mechanisms respectively. It proves that the more the drag force is exerted on the upper and lower sides of the interface the more elongation factor of the centroid is obtained. The average velocities inside the recirculation trapped zone for drainage and imbibition mechanisms were 12.36 µm/s and 10 µm/s respectively. This shows that the presence of two stagnant zones during the imbibition mechanism leads to less intensified recirculation. In other words, less interface is exposed to the high velocity invading phase during imbibition displacement.

For a better understanding of recirculation behavior near the solid surface, drag-force distribution along the right side of the solid circle (close to the interface) inside the left-hand cavity (Figure 10) was analyzed and is shown in Figure 13. The remarkable point is that immediately before and after the center of circulations, two stagnant zones with minimized drag force exerted on the solid circle exist for both imbibition and drainage mechanisms. The dragforce values were higher for the lower part of the semicircle as the radius of the channel around the circle was lower leading to higher velocities. For the drainage mechanism, the drag force suddenly increased at the central part of the semicircle, where the thickness of the trapped layer was significantly decreased. The drag force applied on the semicircle, in the case of the drainage mechanism, was higher as the average velocities inside the trapped zones were 22.46 and 13.55 µm/s, for drainage and imbibition displacements respectively. Analysis of drag-force distribution along the interface for bottom dead-end depicted in Figure 14 shows similar behavior to the other two cases. Although the recirculation of the invading phase in the vicinity of the trapped phase provides enough shear stress to recirculate the trapped phase, this was not sufficient to cause intensified circulation. Analysis of average velocities inside the trapped phase reveals that the velocities inside the trapped phase at the bottom deadend were two orders of magnitude less than the cavities. The elongation factors for centroids of rotation were 0.66 and 1.44 for imbibition and drainage respectively. This is consistent with the results for the drainage mechanism in right-hand cavity. The more drag force exerted on two elongated tails of trapped phase on the solid surface (on which the trapped phase is adhered) results in higher elongation factor of the centroid. The detailed average velocities of trapped phases are shown in Table 3, where the effects of displacement mechanism and fluid configuration on average velocities are obvious. The higher average velocity inside the trapped phase is seen in the case of drainage compared to imbibition.

The average velocity was at a maximum for the left-side cavity where the narrowness of the throat resulted in high momentum transfer from the invading phase to the trapped one. The minimum average velocity inside the trapped phase was at the bottom dead-end, which was behind the mean flow direction. Consequently, the least momentum was transferred to the trapped phase. 4 for various cases. This parameter can give a good sense of the dispersion of the drag force through the trapped phase. Although the values of Vtrap/Aint are slightly higher for the drainage mechanism, the more intensified recirculation can be seen in these cases. This is because of greater accessibility to the interfacial area during drainage as no stagnant zones exist near the interface. The point that draws attention to itself is that the Vtrap/Aint was lowest in the bottom dead-end alongside the minimum average velocity.

As the bottom dead-end was behind the mean flow direction, the drag force exerted on the interface was low. This decreased the average velocity inside the trapped phase in bottom deadend in comparison to two other locations in pore assembly even in the case of low Vtrap/Aint.

Effect of Fluid Dynamic Properties on Recirculation Behavior

Although experimental and numerical studies have analyzed the role of capillary forces and viscosity ratio on the mobilization of the trapped phase as well as the recirculation intensity inside the trapped ganglia [START_REF] Kashid | Liquid-liquid slug flow in a capillary: an alternative to suspended drop or film contactors[END_REF][START_REF] Kashid | Internal circulation within the liquid slugs of a liquid-liquid slug-flow capillary microreactor[END_REF][START_REF] Liu | Micro-PIV investigation of the internal flow transitions inside droplets traveling in a rectangular microchannel[END_REF][START_REF] Ma | On the flow topology inside droplets moving in rectangular microchannels[END_REF]. However, the effects of parameters incorporating in capillary number of invading and trapped phases ) from 1/100 to 100 were considered during an imbibition displacement through the pore-doublet model where a non-wetting phase was patched only on the right-hand cavity. The constant injection velocity was 1×10 -4 m/s and surface tension of 0.05 N/m. As the vertical movements of centers of rotations for different viscosity ratios are negligible, the changes in the horizontal location of centroids to the basis case (M=1) are shown in Figure 15. By deviating the viscosity ratio from unity to the higher or lower ones, the centroids shift to the right. The more the invading viscosity, the more drag force is exerted on the interface, pushing the interface to the right side. On the other hand, as the viscosity of the trapped phase increases, the more energy is dissipated through the trapped phase [START_REF] Pilotti | Energy dissipation and permeability in porous media[END_REF]. Consequently, although the interface moves to the left (due to the change in drag-force balance) the dissipation results in the movement of centroid to the right. In Figure 16 the change in the average velocity of the trapped phase versus the viscosity ratio is shown. It illustrates that by increasing the viscosity ratio, the average velocity inside the trapping phase logarithmically increases. This logarithmic behavior can be attributed to the coupled effects of changes in transmitted shear stress, dissipation of energy inside the trapped phase and movement of centroid.

Figure 15 Changes in the location of centroids for different viscosity ratio

Interfacial Tension (IFT)

As the interfacial tension affects the shape of the liquid/liquid interface [START_REF] Brackbill | Dynamic immiscible displacement mechanisms in pore doublets: theory versus experiment[END_REF], it can directly influence the recirculation of the trapping phase. Therefore, the displacement process from low values for surfactant/hydrocarbon systems to high values for water/hydrocarbon systems is simulated in ranges of 1 to 50 mN/m during an imbibition displacement through the pore-doublet model where the non-wetting phase is patched only on the right-hand cavity. The constant injection velocity is equal to 1×10 -4 m/s and viscosities are 0.01 and 0.001 Pa.s for invading and trapped phase respectively. As the interfacial tension decreases, the curvature of interface changes from a concave shape to a flatter form; consequently, the centroid of recirculation shifts toward the right side of the pore assembly that can be seen in Figure 17. As the interfacial tension increases, the capillary forces are strengthened; as a result, resistance to flow appears which leads to weakened circulation.

Therefore, as it is shown in Figure 18 an overall logarithmic decrease is seen in the average velocity trend as the interfacial tension increases. Changes in interfacial tension reshape the

Injection Velocity

To examine the effect of injection velocities on the velocity field through the trapping phase, various injection velocities from 1×10 -3 to 3×10 -6 with a stepwise increase for imbibition mechanism were considered. Surface tension is equal to 0.05 N/m and viscosities are 0.01 and 0.001 Pa.s for invading and trapped phase respectively. The results illustrated while the increase in velocity field is linearly proportional to the injection velocity can be seen, no change in location of the centroid is observable. Figure 19 depicts the log-log plot of average velocity inside the trapped phase versus the injected velocity. Linear behavior between the average velocity inside the trapped phase and injection velocity can be seen. ) as a function of capillary number of invading phase for various viscosity ratios is evaluated. As it is shown in Figure 20, a linear relationship between the capillary numbers of trapped and invading phases was found, for all viscosity ratios. For the viscosity ratios more than unity (M>1), increase in viscosity of invading phase results in less capillary number in trapped phase as the less drag force is transmitted to the interface. On the other hand, for the viscosity ratios less than unity (M<1), the amount of drag force transmitted to the fluid-fluid interface is constant and increase in viscosity of trapped phase results in a minor increase in capillary number of trapped phase. This is attributed to the fact that the effect of increase in viscosity of trapped phase on the capillary number is cancelled out with intensified energy dissipation inside the trapped phase.

Figure 20 Capillary number of trapped phase as a function of the capillary number of invading phase for various viscosity ratios from 1/00 to 100 (the solid lines are the fitting lines)

Summary and Conclusion

The numerical simulation of the immiscible displacement process through a pore-doublet model was performed in such a way to set trapped-phase patches in specified locations of the pore assembly. To capture the pore-scale processes at low capillary number we used a filtered surface formulation of direct pore-scale modeling for drainage and imbibition mechanisms.

 The results showed, in the case of ranges capillary number of invading phase analyzed in this study (in order of 1×10 -2 to 1×10 -7 ), that a recirculation phenomenon exists in the  Drag-force analysis on the liquid-liquid interface shows that the maximum force was exerted near the center of the interface during imbibition, whereas during drainage the greater force was applied on two elongated interface tails on the solid surface. The elongation factor of the centroid is a function of displacement mechanism (Drag-force distribution along the interface), elongated parallel to the interface during drainage, perpendicularly during imbibition.

 The existence of a solid surface accentuated the recirculation process so that one or more centroids can be created, depending on the displacement mechanism. Analysis of dragforce distribution on the solid surface showed that near the centroids adjacent to the solid surface, stagnant zones were created.

 In the case of the lower ratio of the trapped phase's radius to the cavity depth, two simultaneous recirculation zones inside the adjacent invading and trapped phases were formed.

 Fluid dynamic properties play a pivotal role in recirculation behavior: for all viscosity ratios a linear relationship between the capillary numbers of trapped and invading phases exists.

 The combined effects of changes in the drag-force distribution along with the interface of two-phases as well movement of centroids cause logarithmic behavior of average velocity inside the trapped phase, while this parameter shows a linear trend with changes in injection velocity as the centroid is fixed for various injection velocity.

 These results can be considered as a cornerstone for more complicated cases like CO2 sequestration or reservoir in situ combustion, where the other parameters can be combined with momentum transfer such as mass and heat transfer. Our quantified recirculation results provide great insights into the mobilization process in porous media where forces such as viscosity, inertia, and capillary forces are working simultaneously.

  Figure 1-General Sketch of the problem

  is attributed to the inertial force, which is responsible for the acceleration of fluid and momentum transfer caused by advection. 𝛁. 𝑇 = 𝛁. (𝜇𝛁𝑼) + 𝛁𝑼. 𝛁𝜇 is the viscous force, 𝜌𝒈 is gravity force and 𝒇 𝑐 is the capillary force. 𝑃 𝑑 is dynamic pressure and defined as: 𝑃 𝑑 = 𝑃 -𝑃 𝑐 , where the 𝑃 𝑐 is capillary pressure term. Dynamic pressure is acquired from the mass balance equation in combination with the momentum balance equation, considering incompressible flow: 𝜵. 𝑼 = 0 2

  [START_REF] Heshmati | Interfacial boundary conditions and residual trapping: A porescale investigation of the effects of wetting phase flow rate and viscosity using microparticle image velocimetry[END_REF] study on velocity fields and shear stress at fluid-fluid interfaces in a controlled environment. Figure2shows the 2D geometrical model used in this study. The width of the inlet and outlet of the pore assembly was 120 µm and the model was 410 µm long. This specific model enabled us to examine recirculation inside the trapped phases in various geometrical structures where dead-ends and cavities were designed to resemble the likeliest zones in porous media for the trapping process. The left and right sides of pore-doublet models are cavities perpendicular to the main flow, while in the left side a circular solid provides an opportunity to analyze the rotational flow in the vicinity of solid surfaces. The bottom deadend is behind the main flow direction, allowing investigation of the recirculation in a zone where the invading phase velocities are relatively low.

Figure 2

 2 Figure2Schematic of the pore-doublet model inspired from[START_REF] Heshmati | Interfacial boundary conditions and residual trapping: A porescale investigation of the effects of wetting phase flow rate and viscosity using microparticle image velocimetry[END_REF]

Figure 3

 3 Figure 3 Demonstration of mesh inside (a) the pore assembly, and (b) the left-hand cavity

Figure 4 Figure 5

 45 Figure 4 Comparison of phase distribution and velocity profile for a trapped phase inside the

Figure 6 Figure 7 Figure 9 Figure 10

 67910 Figure 6 Velocity profiles and fluid distributions through the pore-doublet model during, Ca≈3.3×10 -5 , M=10: (a) Imbibition mechanism, and (b) Drainage mechanism (The blue and

Figure 11

 11 Figure 11 Velocity field through the recirculation zone inside the trapped phase at the bottom

Figure 12

 12 Figure 12 Drag-force distribution along the interface per unit of depth (µN/m) at the right-

Figure 13

 13 Figure 13 Drag-force distribution along the right side of the solid circle surface per unit of

Figure 14

 14 Figure 14 Drag-force distribution along the interface per unit of depth (µN/m) at bottom dead-

(

  viscosity ratio on the movement of the center of rotation and average velocity of the recirculating trapped zone are not well discussed in the available literature. Here we analyze recirculation behavior under various dynamic conditions.4.4.1. Viscosity Ratio[START_REF] Li | Numerical simulation of steady flow past a liquid sphere immersed in simple shear flow at low and moderate Re[END_REF] investigated the effect of viscosity ratio at a constant Reynolds number on the velocity field inside a spherical droplet exposed to shear flow. They observed that depending on the viscosity ratio a change in location of the source points (centroids of recirculation) can be seen. In this case, we investigated the effect of viscosity ratio on the recirculation performance at a low capillary number. To understand the role of viscosity ratio, diverse viscosity ratio ranges for various cases including water/polymer solution-heavy oil displacement in porous media (𝑀 = 𝜇 𝑖 𝜇 𝑡

  invading and trapped phases; consequently, the distribution of drag force along the interface is varied. On the other hand, changes in the location of the centroids in conjugation with changes in drag force distribution cause the logarithmic (non-linear) behavior of average velocity inside the trapped phase.

Figure 17 Figure 18

 1718 Figure 17 Changes in the location of centroids for different interfacial tension

Figure 19

 19 Figure 19 Average velocity of trapped phase versus injection velocity 4.4.4. The Role of Capillary Number

  while various behaviors were observed. Zero-velocity zones were observed at the centroids of the recirculation trapped phase.  In the imbibition mechanism, the convexity of the trapped phase created two stagnant regions in the vicinity of the interface inside the invading fluid. In the drainage mechanism, the concavity of the trapped phase wiped out the stagnant zones inside the invading fluid adjacent to the fluid-fluid interface.

Table 1

 1 Fluid properties

	Fluid type	Viscosity (Pa.s)	Density (kg/m 3 )	Surface tension
				(N/m)
	Invading phase	0.01	1000
				0.03
	Trapped phase	0.001	1000

Table 2

 2 Comparison of velocities near the interface inside the invading and trapped phases in experimental and modelling approaches

		Velocity inside	Velocity inside non-
		wetting invading	wetting trapped
				Velocity ratio
	Approach	phase (10 µm far	phase (10 µm far
				(invading to trapped)
		from the interface)	from the interface)
		(m/s)	(m/s)
	Experimental results	110 (µm/s)	21 (µm/s)	5.2
	(Heshmati and Piri,
	2018)
	Modelling results	99 (µm/s)	23 (µm/s)	4.3
	(this study)

Table 3

 3 Average velocity inside the trapped phase for different fluid configurations and Another factor that is of great importance is the volume of the trapped phase normalized by interfacial area, Vtrap/Aint, listed in Table

		various displacement mechanisms
	Displacement	Average Velocity inside the recirculating trapped zone (µm/s)
	Mechanism	Right-side Cavity	Left-side Cavity	Bottom dead-end
	Drainage	12.4	22.5	0.7
	Imbibition	10.0	13.5	0.1

Table 4

 4 Volume of the trapped phase normalized by interfacial area for different fluid

	configurations and various displacement mechanisms

Acknowledgments

We gratefully acknowledge the financial support provided to the PIVOTS project by the "Région Centre -Val de Loire" and the European Regional Development Fund.

Appendix A

𝜌 and 𝜇, respectively density and viscosity of the VOF mixture, are calculated using the indicator function 𝛼 ∈ [0, 1] representing the volume fraction of one of the phases in a controlled volume:

𝜌 = 𝛼𝜌 𝑤 + (1 -𝛼)𝜌 𝑛𝑤

A-1

where the indices w and nw denote wetting and non-wetting fluids respectively.

Inaccurate representation of capillary forces results in smeared color function. [START_REF] Raeini | Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method[END_REF] suggested a sharpening color function to solve this problem. Using Sharp Surface Force modeling (SSF) approach, capillary body forces based on Equation 6 are given by,

where the smoothed color function is given by,

where the indices center and face denote the cell center and face center respectively. Also, the modified sharpened color function is given by [START_REF] Raeini | Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method[END_REF],

where 𝐶 𝑠ℎ𝑝 ∈ [0,0.5) denotes the heuristic sharpening coefficient. While 𝐶 𝑠ℎ𝑝 = 0 results in original CSF formulation, as 𝐶 𝑠ℎ𝑝 approaches 0.5 becomes sharper.

Although the SSF formulation can appropriately eliminate the non-physical velocities in the absence of complex solid boundaries, the parasitic currents appear close to curved or edged solid boundaries. To address this issue, capillary-induced flows that results in currents parallel to fluid interfaces are modified as:

The capillary forces parallel to the interface starts from zero as initial condition and is updated as follows :

𝒇 𝑐,𝑓𝑖𝑙𝑡∥ = 𝛿 1 (𝒇 represent the capillary term parallel to the interface.

Appendix B

The figures below show how the velocity is distributed inside the trapped phases for various locations of pore assembly. Also, it shows that the velocity inside the trapped phase is noticeably less than the adjacent invading phase. Closing to the wall solid surface the velocity decreases. The direction of rotation inside the trapped phase is a function of fluid configuration within the pore assembly as well as the direction of flow.