Adaptive mesh refinement in locally conservative level set methods for multiphase fluid displacements in porous media

Multiphase flow in porous media often occurs with the formation and coalescence of fluid ganglia. Accurate predictions of such mechanisms in complex pore geometries require simulation models with local mass conservation and with the option to improve resolution in areas of interest. In this work, we incorporate patch-based, structured adaptive mesh refinement capabilities into a method for local volume conservation that describes the behaviour of disconnected fluid ganglia during level set simulations of capillary-controlled displacement in porous media. We validate the model against analytical solutions for three-phase fluid configurations in idealized pores containing gas, oil, and water, by modelling the intermediate-wet oil layers as separate domains with their volumes preserved. Both the pressures and volumes of disconnected ganglia converge to analytical values with increased refinement levels of the adaptive mesh. Favourable results from strong and weak scaling tests emphasize that the number of patches per processor and the total number of patches are important parameters for efficient parallel simulations with adaptive mesh refinement. Simulations of two-phase imbibition and three-phase gas invasion on segmented 3D images of water-wet sandstone show that adaptive mesh refinement has the highest impact on three-phase displacements, especially concerning the behaviour of the conserved, intermediate-wet phase.


Introduction
Multiphase flow in porous media occurs in a wide range of engineered applications in the subsurface, including enhanced oil recovery, permanent CO 2 storage, temporary gas storage, and remediation of non-aqueous pollutants.For all these applications it is important to determine the extent of bypassed or trapped fluids after invasion of another fluid.This Computational fluid dynamics (CFD) techniques for modelling multiphase flow at the pore scale include lattice-Boltzmann methods [29,30,35,49,57,59,67,71], smoothed particle hydrodynamics methods [32,65], and phase-field methods [31,60,72].More conventional CFD methods typically combine the Navier-Stokes equations with a technique for interface capturing, like the volume-of-fluid (VoF) method [23,48,56] or the level set (LS) method [44,55].These interface-tracking methods are implicit, so that they can handle topological changes (like ganglia splitting and merging) in a natural manner.The VoF method describes the phase locations with an indicator function, and its motion is governed by a conservation law [23].Thus, the VoF method has intrinsic volume/mass conservative properties, but it requires advanced numerical methods to handle the sharp changes of the indicator function at interfaces and for calculating interface properties like curvature.The LS method is based on advection of a smooth LS function typically given as the signed distance to the interface.As such, the LS method is not conservative, but interface properties like curvature is easily calculated directly from the LS function.Several techniques has been developed to improve conservation of LS methods, including the implementation of higher order numerical schemes [40], particle LS methods [9], combined LS and VoF methods [61], the conservative LS method which adopts a modified LS function [41], as well as various interface-correction methods [13,37].
Using the above methods for direct simulation of multiphase flow on segmented images of the pore structure, while accounting for capillary, viscous, and gravitational forces, is computationally costly.However, for multiphase flow at the pore scale in reservoir rock, capillary forces normally dominate over viscous forces and gravity [21].Thus, it is adequate to approximate the flow as multiphase displacements controlled by capillary pressure, which makes for more efficient simulations.These displacements are generally modelled as a series of quasi-static equilibrium states for different prescribed capillary pressures [17,69].Such quasi-static simulations generate capillary pressure-saturation curves for drainage and imbibition that can be correlated against experimental data.Methods in this category include pore-network modelling based on invasion percolation [2,43,45,73], poremorphology methods [15,22,39,66], and LS methods [18,27,47].Within interface-capturing methods, we also note that the MBO method [38] (named after its authors) and related variants [51,62], which solves a diffusion equation and updates an indicator function in alternate steps, can be used to simulate curvature-driven motion in multiphase systems.Common to all these approaches is that they need to be supplied with methods for local conservation of mass (or volume, for incompressible fluids) to accurately describe the behaviour of fluid ganglia that detach from the continuous phases.
Jettestuen et al. [28] presented a new method for local volume conservation (LVC) of fluid ganglia within the setting of capillary-controlled displacement, using LS methods for two-phase systems [27] and the multiphase level set method (MLS) for three-phase systems [18], abbreviated as LS-LVC and MLS-LVC methods, respectively.The LS/MLS-LVC methods have the option to enforce conservation on selected phases, and it calculates surface area, volume, and pressure of individual fluid ganglia of the conserved phases during displacement.A volume redistribution algorithm handles volume conservation during splitting and coalescence of isolated domains.Jettestuen et al. [28] applied the LS/MLS-LVC methods to various 2D and 3D geometries, where it showed a quadratic volume convergence with respect to resolution of the uniform grid.The model was used to simulate gas invasion after primary drainage by conserving only the oil phase in one case and both oil and water phases in another.It demonstrated that local volume conservation of both phases predicts a higher residual oil saturation than only oil conservation.However, we expect finer grids or adaptive mesh refinement (AMR) to improve the resolution of this problem, especially with respect to oil layers in three-phase systems and small fluid ganglia.Consequently, the improved model should generate more realistic capillary pressure curves and distributions of trapped fluid ganglia, allowing a more accurate ganglion characterization.Thus, the objective of this work is to extend the LS/MLS-LVC methods such that they can utilize AMR capabilities.
AMR is a computational technique that can help in the dynamic focussing of refinement to resolve local features in the computational domain without incurring time penalties associated with global refinement.It is an effective tool for improving the accuracy of numerical simulations while maintaining reasonable computational efficiency [10].In general, there are two main categories for AMR implementation, the patch-based approach [4,11,14] and the tree-based approach [46,63,70].The patch-based approach requires tagging unresolved regions using special refinement criteria and then grouping these grid cells into patches with a finer mesh [14].On the other hand, a quadtree approach works on the principle of recursive decomposition, where a quadtree stores the data [63].Although the second method is computationally fast and efficiently uses storage memory, its application can require complicated algorithms and interpolations to manage hanging nodes.Jettestuen et al. [28] implemented the LS/MLS-LVC methods for uniform grids using the software framework SAMRAI [3,24,25], which is a C++ library for parallelism and AMR that uses ghost cells and a patch-based approach to provide excellent parallel scalability while handling a large number of processors [14].Ghost cells of a patch are extensions into the neighbouring patches where the values on the grid are generated through coarsening and refinement of values on the neighbouring patches.
Coarsening operations transfer data from finer grids to coarser grids, and refining operations transfer data from coarser grids to finer grids.
In this work, we extend the work of [28] by including AMR capabilities for the LS/MLS-LVC methods and integrating it in the parallel framework of SAMRAI.While SAMRAI provides AMR capabilities for real numbers (e.g., to describe scalar fields), the challenge in our LVC application is handling integers on a map that represents isolated domains.The paper is organized as follows: First, Section 2 briefly describes the LS/MLS-LVC models, while Section 3 presents the methodology and related algorithms for incorporating AMR in the LVC method.In Section 4 we validate the MLS-LVC implementation with AMR against analytical solutions for three-phase configurations in a triangular pore.We discuss the computational efficiency obtained by using AMR rather than a corresponding uniform fine grid in the simulations.Then, we provide weak and strong scaling test results to demonstrate the parallel scalability of the new algorithm.Subsequently, we investigate the impact of AMR on simulations of two-phase imbibition and three-phase gas invasion with disconnected oil ganglia on 3D segmented pore-space images of sandstone.Finally, Section 5 summarizes the results.

Methods
We investigate capillary-controlled multiphase systems of gas (g), oil (o) and water (w) at the pore scale, in which a series of capillary equilibrium states for different imposed capillary pressures describes the fluid displacement.At these equilibrium states, the capillary pressure across the fluid/fluid interfaces in the pore space satisfies the Young-Laplace equation: where P αβ [Pa] is the capillary pressure, σ αβ [Nm -1 ] is the interfacial tension between phases α and β, p α [Pa] is the pressure of phase α, and C αβ [m -1 ] is the interface curvature.The fluid/fluid interface meets the pore/solid interface at a contact angle θ αβ , which is measured through the denser fluid phase β and satisfies Young's equation: where σ sα and σ sβ are the solid/fluid interfacial tensions.By combining Young's equations for all three fluid pairs, we obtain the Bartell-Osterhof equation [e.g.,8], which constrains the interfacial tensions and contact angles in three-phase systems at thermodynamic equilibrium.
In this section, we briefly describe the LS [27] and MLS [18] methods for capillary-controlled displacement, as well as the technique used for incorporating local volume conservation (LVC) in these methods [28].

Level set method for two-phase capillary-controlled displacement
The level set (LS) method is an Eulerian approach to model propagation of interfaces under the influence of a velocity functional [1,44,52].In this method, a level set function φ, normally given by the signed distance function, describes the interfaces as φ = 0, so that φ takes different signs on the inside and the outside of a fluid.The equation of motion for an evolving level set function φ under action of a normal velocity F N and an advection velocity F A is given by Adalsteinsson and Sethian [1]: During evolution with Eq (4), the level set function φ gets distorted and requires repeated reinitialisation to maintain its signed distance nature.The reinitialisation is carried out using the following equation [44]: where S(•) is a smoothed sign function given by and x is the grid spacing.LS methods have previously been used to model capillarycontrolled displacement at the pore scale [27,47].Following Jettestuen et al. [27] the LS function φ describes the interface between fluid phases α and β as the zero contour φ = 0, where φ < 0 in phase α and φ > 0 in phase β.Another LS function ψ is static and describes the pore geometry, so that ψ > 0 in pore space, ψ < 0 in solid, and ψ = 0 on the pore/solid interface.A sharp Heaviside function H (•) separates the velocity term in the pore space, which governs the capillary-controlled displacement, from the velocity term in the solid space, that forms the contact angle θ on the pore walls.The evolution equation (4) takes the following form [27]: Here, At capillary equilibrium (φ t = 0) the velocity terms in equation (7) are zero in vicinity of the interfaces.Then we obtain C = κ (Young-Laplace equation) in the pore space, while in the solid space and on the pore walls, the intersection angle ζ satisfies where n φ and n ψ are unit normal vectors of φ and ψ, respectively.After each calculated equilibrium state using equation ( 7), we can increase or decrease C stepwise to simulate drainage or imbibition, respectively.A more detailed description of the LS method with numerical validation and applications are provided by Jettestuen et al. [27], Helland et al. [17], and Friis et al. [11].

Multiphase level set method for three-phase capillary-controlled displacement
The MLS method describes each fluid phase α by its own unique LS function φ α , so that φ α < 0 inside phase α, φ α > 0 on the outside, and φ α = 0 on the boundary of α to the other phases.While this representation can describe an arbitrary number of phases, we focus here on three-phase displacement where α = g, o, w.For each fluid phase, we also specify a phase pressure p α , surface tension γ α , and a fluid/solid intersection angle ζ α which is measured through phase α and defined by cos ζ α = n φ α • n ψ .The surface tensions and intersection angles satisfy the following equations: where αβ = go, ow, gw, and γ s is an imaginary surface tension for the solid phase that ensures σ αs is always positive.Hence, equation ( 12) is the same as Young's equation (2).Following [18], the intersection angles are taken as which, by equation ( 12), leads to a set of contact angles that fulfills equation (3) for a given set of γ α .
The evolution equations for the fluid LS functions φ α are given by [16,18]: Here, γ and ζ remain constant throughout the simulation.On the other hand, p α is constant, as specified by the boundary conditions, for continuous phases whereas it varies with time for isolated domains.The calculation of isolated domain pressures is detailed in the next section.
As these evolution equations are uncoupled, the three LS functions φ α will develop overlap and void regions around the interfaces during evolution.Consequently, in an interface region of two phases α and β, the intersection of φ α and φ β occurs at a nonzero contour (i.e., φ α = φ β > 0 in case of void region or φ α = φ β < 0 in case of overlap region) rather than at the zero contour (φ α = φ β = 0) which would correspond to a unique representation of the αβ-interface.This problem is solved using the projection method of Lossaso et al. [36] where, in every grid cell, the average of the two smallest φ α is subtracted from all φ α at the end of each time-iteration step.This step moves all intersections of the different LS functions to zero contours, creating a unique delineation of the phases in the computational domain.When the multiphase system has relaxed to a steady state, i.e., (φ α ) t = 0, α = g, o, w, the solution of equation ( 14) (with projection step included) satisfies Young-Laplace equation in the pore space and Young's equation on the pore walls [16,18]: since φ α = −φ β in the vicinity of the αβ interfaces, which implies ∇φ α = −∇φ β and κ α = −κ β .
The MLS method uses non-dimensional parameters for length, surface tension and pressure, L, γ , and p, respectively.These parameters are derived from actual values of length L actual [m], surface tension γ actual [Nm -1 ] and pressure p actual [Pa], using characteristic properties for the physical system: where L * [m], γ * [Nm -1 ], and p * [Pa] are the characteristic parameters.
In both LS and MLS methods, the numerical discretization approximates the normal and advective velocity terms using a weighted essentially non-oscillatory (WENO) scheme with appropriate upwinding techniques, while the curvature terms use central differences.The iteration-time discretization uses a third-order Runge-Kutta method with a time step determined from a standard Courant-Friedrichs-Lewy (CFL) condition [44].In this work, the CFL coefficient is 0.6.Faster temporal discretization methods exist for steady-state problems on structured mesh [33,68], but it remains to investigate their applicability to capillary-controlled displacement on complex pore geometries.The convergence in the MLS method is declared when the difference of φ α between two reinitializations, denoted by m and n, in a small band around the zero contours of the LS functions in pore space is lower than a specified tolerance value: where H is a smoothed Heaviside function and = 1.5× x is the smoothening parameter.In equation ( 18) λ = b × is the width of the convergence band.In this work, we set b = 5 and c = 0.001.
Note that the interface position can slightly shift after a reinitialization step.This effect is most significant in MLS simulations of triple junctions between the three fluids where rounded features in the zero contours could form and result in small overlap/void regions.This impact of reinitialization is most severe when the angle at the triple junction through one of the fluids is small (as determined by the interfacial tensions).In our experience this region is about the size of one to two grid cells, so the impact is small.The projection step [36] in the MLS iteration after the reinitialization removes these small regions.

Local volume conservation method
Jettestuen et al. [28] developed a technique for local volume conservation (LVC) of fluid ganglia in the LS and MLS methods.Here, we describe the LVC method for a selected conserved phase α in the MLS setting.The main idea is to determine a spatially and temporal phase pressure p α (x, t) for use in equation ( 14) that conserves the volumes of disconnected ganglia of phase α during level set evolution.To this end, the method first identifies the disconnected domains of phase α using a grassfire method [47] that is based on Chebyshev distance (i.e., L ∞ -norm) as metric with 8-connectivity (in 2D) and 26-connectivity (in 3D).Note that the Chebyshev distance from the gridcell to these connected blocks is 1 unit.We denote the identified domains at iteration time t as a (t), where a ∈ I α (t), and I α (t) is a numbered list of all ganglia of phase α at time t.Then, we construct extended domains ext a (t) that completely encloses domains a (t), so that ext a (t) ⊃ a (t).The extended domains include all points nearest to the boundary of the corresponding domains as determined by a distance map, using city-block distance (i.e., L 1 -norm) as metric.Note that city-block distance (also referred to as Manhattan distance or L 1 -norm) is the minimum number of grid cells traversed laterally to reach a point from another point.The extended domains are needed to distribute volumes correctly from one iteration time step to another, but they also ensure accuracy in the numerical discretization of equation ( 14) across interfaces.These extended regions partition the computational domain completely: Here, we have distinguished the extended disconnected regions ext a (t) from the extensions ext α,bnd (t) of all regions α,bnd that are connected to inlet or outlet boundaries and referred to as the continuous phase.The continuous phase changes volume as phase α enters or leaves the computational domain.
Then, the phase pressure for a conserved phase α is given by Saye and Sethian [52]: Here, p α,bnd is the continuous-phase pressure assigned to regions connected to inlet or outlet boundaries.If phase α is the invading phase, stepwise increments of p α,bnd after every capillary equilibrium state drives the quasi-static displacement.Further, V (0) a (t) is the target volume of a disconnected region a , while V a (t) and A a (t) are volume and surface area of the region, respectively, calculated at each iteration-time t: where δ is a smoothed Delta function.Calculation of global saturation uses equation (21) with integral over the entire computational domain .These smoothed Heaviside and Delta functions are calculated as: and, A conservative volume redistribution technique assigns volumes to disconnected domains at the next time step based on the volume of domains at the previous time step.The method accounts for potential ganglia splitting and coalescence.Assume we have identified the sets of domains and extended domains at iteration time t The target volume assigned to domain b (t 2 ) is given as where | • | denotes the volume measured in number of voxels.Equation ( 26) uses the relative volume of overlap between regions at t 1 and t 2 as a weight for the volume redistribution.The set of extended regions cover the entire computational domain, see equation (19), which ensures volume conservation during the redistribution.Further, restrictions on time step in the MLS method guarantees small interface changes within an iteration step, so that target domain volumes are redistributed correctly.We remark that continuous-phase domains ext α,bnd (t) are not part of the redistribution equations as their pressure is prescribed.Their role is merely to limit the extensions of nearby isolated domains.When new regions form by detachment from the connected phase, their target volumes V (0) a (t) are calculated from equation (21).Thus, equation (20) implies that the pressures of such regions will be zero in the first time iteration step.However, their temporary volumes V a (t), also calculated by equation (21), will generally differ in subsequent iterations, resulting in nonzero domain pressure.
The phase-by-phase representation of the MLS method allows for applying LVC to one, two or three phases in a three-phase system.The LS method, which instead uses one LS function to describe interfaces in a two-phase system, will also require a similar phase-by-phase representation to account for LVC of both phases.However, conservation of one phase follows the same approach as described above.The main difference is that the LS method constructs a spatial and temporal interface curvature C(x, t) that replaces pressure (equation ( 20)): where ϒ = 1 if φ < 0 represents the conserved phase, and ϒ = −1 if instead φ > 0 represents the conserved phase.

LVC methodology for AMR and parallelism
Here we introduce the methodology and algorithms for parallel AMR implementation of MLS-LVC method.Apart from the differences described in the previous section, the AMR implementation for LS-LVC method will be similar.An important part of the LVC algorithm presented by Jettestuen et al. [28] is the identification of isolated fluid domains.The software implementation uses a patch-based approach, so the identification must first be performed on each patch.However, since a domain usually crosses (often several) patch borders, unique global labelling (using integers) is necessary.As discussed in [28], the corresponding computational algorithm requires considerable communication in its parallel version, which is implemented using Message Passing Interface (MPI).The LVC algorithm constructs extended domains by extending the global labels outward into the computational grid from the interface of each isolated domain.
In the parallel setting, patch boundaries limit these extensions, so that equation (19) does not generally hold.Thus, additional extension procedures are carried out across patch boundaries to construct sufficiently large extended domains.As mentioned in Section 2.3, the role of the extended domains are as follows: (i) to conserve domain volumes during evolution, (ii) to redistribute volume during domain splitting and coalescence, (iii) to calculate volume and surface area of the domains, and (iv) to distribute the calculated phase pressure on the computational grid (see also sections 3 to 5 in [28] for details on LVC application for uniform grids).
A major challenge in applying local volume conservation on a uniform grid is creating a computationally efficient way of identifying isolated domains and populating appropriate domain pressures on the geometry.Furthermore, these pressures should be suitably extended from the interfaces so that the size of the numerical stencil does not affect the level set evolution.We handle this requirement by linking pressures to the extended regions and using city-block distances for modifying extended regions to suitable distances from the interfaces across patches.However, the case for AMR grids is more complex than for uniform grids.Firstly, AMR grids require additional attention in handling integers during data communication between patches of different sizes.Secondly, the city-block distance between patches with different grid spacings will not correspond to each other.Finally, AMR application requires regular regridding (changing the grid structure) to maintain accuracy and be computationally efficient.Consequently, additional efforts are needed to ensure that various variables involved keep the same values before and after regridding.
The MLS-LVC model is implemented in C++ programming language using the SAMRAI software framework.SAMRAI allows flexibility in cell refinement criteria and frequency of regridding.The gridding and regridding process is accomplished in a linear order time [14].We have utilised the framework to generate patches with finer mesh and communicate relevant values between patches through coarsening and refining variables in the ghost cells.In SAM-RAI, adaptive grids are a layered structure where finer grids are successively placed on top of the coarsest grid (see Fig. 1).The calculations are done on the topmost layer or mesh (i.e., the finest grid) for each overlapping location.When needed, the values are transferred to other meshes by coarsening or refining operations.
In the following subsections, we provide the LVC algorithm for AMR grids followed by a discussion on the necessary modifications made to the original (uniform grid) algorithm for AMR compatibility.

MLS-LVC algorithm with AMR
The LVC algorithm for AMR grids is a slight modification of the parallel algorithm for uniform grids from Jettestuen et al. [28](Section 5.1).The improvements reduce the memory space requirement for uniform grids while adding the capability to handle AMR grids.The steps in the modified algorithm are as follows, assuming α is a conserved phase: 1. Initialize the system at initial iterative time t 0 : Loop over iteration time (from t n to t n+1 ): 2. Set pressure of conserved phase, p α (x, t n ).For continuousphase domains, assign phase pressure p α,bnd to ext α,bnd (t n ).For isolated domains, set pressures in ext a (t n ), a ∈ I α (t n ), using equation ( 20) (enforcing volume conservation) as follows: (a) On each patch: Perform computations according to equations ( 21) and ( 22).(b) On each processor: Aggregate the volume and area information and communicate it to the root processor.(c) On the master processor: Calculate pressure of isolated domains using equation ( 20) and communicate relevant information back to all processors and patches.
3. Iterate the multiphase system from t n to t n+1 with equation ( 14) using the phase pressure from the previous step and then enforcing the projection step [36].This is done in parallel using the patch-based parallelism in SAM- RAI.Here, we also periodically reinitialise level sets and check for convergence.4. Repartition the system at time step t n+1 , i.e., construct isolated domains { a : φ α ≤ 0 and a − b ∞ > 1, ∀ a = b , a and b ∈ I α (t n+1 )} in a similar manner as in step 1 (a) above.Then, extend the isolated domains to obtain the set { ext a (t n+1 )} in a similar manner as in step 1 (b) and (c) above.5. Redistribute volumes according to equations ( 25) and (26) while excluding continuous phase domains ( ext α,bnd ). 6. Identify new regions detached from the continuous phase and calculate their region target volumes based on equation (21).7. Regrid at regular intervals {n mod Q = 0 : Q is the regridding interval} (Algorithm 4, Section 3.2).8. Go to step 2 with n ←− n + 1.

Considerations for adaptive mesh refinement
Transforming the LVC algorithm for compatibility with an adaptive grid necessitates reconciling variables between meshes with different grid spacings.Compared to integers, handling real numbers is easy; the values from a grid of one size can be easily transferred onto a grid of another size through interpolation, extrapolation, or averaging.SAMRAI provides operators that carry out these operations for real numbers.However, averaging integers can produce undesired results for our purpose.For example, if four adjacent fine grid cells have values 1, 2, 4 and 5 in the extended domain map, then the coarser grid cell below them will be wrongly filled with 3 after averaging.Thus, we define our own operators for coarsening and refining of the integer data type (see Fig. 3).
Refining is easy since it creates copies of the coarser value on the finer grid (see Fig. 3 (a)).Coarsening is more challenging since the finer grid cells over a coarser grid cell can either have equal values or different values.So, we create an operator that copies the value on to the coarser grid cell if all the finer grid cells over it have the same value (see Fig. 3  (b)), otherwise it returns the minimum possible C++ integer value, INT_MIN (see Fig. 3 (c)).This unique value on the coarser grid allows us to check if the overlaying finer grid cells have the same values or not.We use this operator on patch variables M and E that store the isolated domains map and the extended domains map, respectively.This operator helps in preserving most of the original values on the coarsened regions and the rest can be recovered by extending the map values on the new grid using Algorithm 2 and 3. We also create another coarsening operator that always returns INT_MIN.This second operator is used for dynamic identification of refined domains on the coarser patch throughout the simulation.In this way the various computations (e.g., volume computations) can be performed with good parallel efficiency.We use this second operator on patch variable P that stores the processor rank.Since Open MPI library provides only positive integers as processor identification number, INT_MIN acts as a distinct value that will not be provided by the MPI library.
In our model, integers are used for identifying isolated domains a and their extended domains ext a , a ∈ I α (t).Algorithm 1 describes the steps followed in identifying isolated domains of a particular phase on the computational grid in the modified algorithm.We identify isolated domains on a given patch using the grassfire algorithm with φ α as the input.The grassfire algorithm (referred to as grassfire transform in image processing) can be described as 'setting fire' to the zero level set contour.It finds the closest region to each grid cell in the patch.The algorithm requires alternate sweeps of the patch from any two diagonally opposite corners while measuring backward distances from the area already swept.In this patch-level identification, we use INT_MAX (the maximum integer limit in C++) to represent domains connected to the computational grid boundary.We label the isolated domains found on each processor sequentially (starting from 1) but the uniquely labelled isolated domains on a processor are not necessarily disconnected from each other.Thus, we use the processor number and the local domain number to identify each domain uniquely.

Algorithm 1 Identify domains.
Require: Patch variables: Phase level set value (φ α ), processor rank (P), domain map (M), and number of grid cells in z-direction (n z ), y-direction (n y ), and x-direction (n x ) Ensure: Steps 6 to 18 of Algorithm 1 check if two identified isolated domains are adjacent to each other.Since step 4 prevents such an identification inside the patch, it is sufficient to look at patch boundary grid cells and the ghost cell adjacent to them to find such pairs.Steps 9 to 17 of Algorithm 1 ensure that for each neighbouring pair spanning adjacent patches, one patch identifies the pair and the other ignores it.We analyzed and compared the computational efficiency of sending the complete initial dataset of identified neighbours to the root processor vs. sending a reduced dataset to the root processor after identifying unique isolated domains on each proces-

Algorithm 2 Define extended domains
Require: Patch variables: Isolated domain label (M), extended domain label (E), distance (D), ratio of patch grid spacing w.r.t.finest patch grid spacing (R), and number of grid cells in z-direction (n z ), ydirection (n y ), and x-direction Update E i, j,k to E of grid cell that modified D i, j,k above 11: end if 12: end for 13: Update E i, j,k to E of grid cell that modified D i, j,k above 17: end if 18: end for 19: In case of AMR simulations, repeat steps 4 to 12 20: end if sor.Our analysis shows that sending the original dataset to the root processor provides faster results.This study is detailed in Appendix A. On the root processor, steps 21 to 25 in Algorithm 1 provide globally unique labels to isolated domains found on the patch (for details refer [28], Appendix A). Figure 2 (top row) depicts an example of isolated domain identification with different levels of refinement.

Algorithm 3 Modify extended domains across patches
Require: Patch variables: Processor rank (P), extended domain label (E), distance (D), ratio of patch grid spacing w.r.t.finest patch grid spacing (R), modification distance threshold (d c , in city-blocks), and number of grid cells in z-direction (n z ), y-direction (n y ) and x-direction (n x ) 1: if region with values from refined areas exists (identified by P = INT_MIN) then 2: end for 8: α ), temporary volume (V α ), and temporary area (A α ) to store V (0) a , V a , and A a (where a ∈ M), respectively for use in equation (20).Ensure: New variables are sized according to the new grid geometry 1: Reinitialize the phase level set values to signed distance function 2: Tag cells for refinement and regrid 3: Refine and coarsen ghost cell data 4: Update the extended domains using Algorithm 2 5: Modify the extended map using Algorithm 3 6: Create new variables for volume and area, V (0) * α , V * α , and A * α , sized according to new grid geometry 7: Fill the new variables (V (0) * α , V * α , and A * α ) using old variables (V (0) α , V α , and A α ) such that E new = E old 8: Free memory used by variables (V (0)  α , V α , and A α ) based on old grid geometry and use the new variables for further pressure calculations We generate the set of extended domains ext a using Algorithm 2 on each patch.On uniform grids, we accomplish it by carrying out two sweeps of the patch for each grassfire run (steps 4 to 18).But for AMR grids, the extended regions are generated by running three sweeps of the patch.With AMR, some parts on coarser patches have INT_MIN values since finer patches cover them.Generally, in our model, the finer patches are not distributed in such a complicated manner that multiple grassfire run sweeps are required to distribute the pressure values effectively.In such scenarios, the additional sweep of the patch (step 19, which is a rerun of the first sweep) provides the correct minimum distance needed in the grassfire algorithm.Figure 2 (middle row) shows the result after extension with different levels of refinement.
These extended regions on the patch are further extended into neighbouring patches by using distances from interfaces and patch neighbourhood information (see Algorithm 3).As mentioned earlier, the city-block distances are measurements where one unit distance is equal to one grid cell.However, the city-block distance on patches with different grid spacing will represent different actual distances.To ensure that this difference does not lead to insufficient domain extensions across borders of patches with different grid spacings, we use the ratio of patch grid spacing to make the distances proportionate such that the finest grid cell equals one unit of distance.For example, suppose the coarse grid spacing is two times bigger than the finest grid spacing.In that case, one fine-patch grid cell equals one unit of distance, and one coarse-patch grid cell equals two units of distance.Another point of importance is the order in which these extensions are carried out from one patch to another.For uniform grids, the order of extension is inconsequential.However, on AMR grids, the order is important to ensure that the level set velocities near zero contour are independent of the patch structure.So, we first carry out extensions from finer patches to coarser patches below them (steps 1 to 8 of Algorithm 3).This first run ensures that if an interface lies only on the finest patch, its effect extends into all the neighbouring coarse patches.Then, we carry out a second extension (as for uniform grid algorithm in Jettestuen et al. [28]) on all the patches (steps 9 to 14).Since the minimum required depth of modification depends on the numerical stencil size, an optimal way to implement this algorithm is to find the necessary patch boundary grid cells.In steps 1 to 8, this boundary is the boundary of a finer patch lying inside the coarse patch, and the distance of the ghost cell underlying the neighbouring finer patch from the closest domain on that finer patch is D f iner−ghostcell .While for steps 9 to 14, it is the boundary with a patch on the same or lower level, and the distance of the ghost cell from the closest isolated domain on that neighbouring patch is D ghostcell .Once the requisite boundary grid cell has been found, we can check whether the extended domain value on the grid cells within the required distance from this grid cell needs to be modified.Figure 2 (bottom row) shows the result after modification of extension maps with different levels of refinement.
In our model, we store the original volume, temporary volume, and temporary area (required in equation ( 20)) in 1D arrays based on the patch grid structure.We number the patches on each processor sequentially and store these variables for extended domain labels present on that patch.For example, if patch number 1 on processor 0 has extended domain labels 2 and 3, then our array starts with the original volume, temporary volume, and temporary area for domains 2 and 3, respectively.If patch number 2 on processor 0, has extended domain labels 1, we append the original volume, temporary volume, and temporary area for domain 1 to these arrays.This array shape helps in limiting the communication between root processors and other processors since we can use M P I _Gatherv and M P I _Scatterv to share the arrays efficiently.We update these variables after each regridding operation using Algorithm 4.During regridding, we create similar variables based on the newer grid configuration.We refine and coarsen all patch variables using integer operators (described previously) and real operators (provided by SAMRAI) before any other operation.Then, we extend the domains using Algorithms 2 and 3. We utilize the extended domain labels to link the variables based on the previous grid structure to those based on the new grid structure.Then, we erase the older variables from memory space and use the newer ones in their place for further calculations.

Results and discussion
Here, we present simulation results with the LS/MLS-LVC models using the new AMR capabilities.We start by validating the model against analytical results for three-phase configurations in a 2D triangular pore and a 3D cylindrical pore, followed by analysing the parallel performance of the AMR implementation through scaling tests.Then, we apply the LS/MLS-LVC models on segmented micro-CT images of a sandstone to explore the significance of AMR in simulations of two-phase and three-phase capillary displacements with disconnected fluid ganglia.
In this section, the phrase '1-level' or 'AMR-1' describes simulations using AMR with one level of refinement, while '2-levels' or 'AMR-2' describes simulations using AMR with two levels of refinement.The phrase 'LVC-1' and 'LVC-2' refer to local volume conservation of one and two phases, respectively.For visualization purposes, we use the software VisIt developed by LLNL [6].To visualize the fluid configuration inside a pore geometry, we confine the zero contours of φ α to the pore space using φ α = max(φ α , −ψ), α = g, o, w.

Validation against three-phase fluid configurations in a 2D triangular pore
We consider three-phase fluid configurations in a 2D triangular pore where water is the wetting phase residing in pore corners, oil is present as intermediate-wet layers, and gas is the non-wetting phase occupying the central portion of the triangle.This configuration is similar to the one theoretically studied by Hui and Blunt [26].We will assume the oil layers are disconnected from each other with different pressures and their individual volumes (areas in 2D) conserved.This situation is different from our previous investigations on triangular pores using MLS method where global oil volume conservation led to a common pressure for all oil layers [18].For stepwise increments of gas pressure at constant water pressure we simulate the displacement of oil layers toward pore corners, applying LVC on the oil phase only.Our objective is to show that higher levels of AMR improve the accuracy in simulations based on a comparison with analytical solutions and that AMR simulations provide similar results as simulations on a corresponding uniform, fine grid.
Geometrically, an oil layer can exist in a corner if the contact angles and interface curvatures, θ αβ and C αβ , αβ = go, ow, satisfy the relations [26]: and where k is the half-angle of corner k.Here we assume that an oil layer ceases to exist when the gas-oil and oil-water interfaces form a triple junction on the pore walls or at the corner bisector.The volume V o of an oil layer in corner k, sandwiched by corner water and bulk gas, is [18] V where From equations ( 30) and ( 31) we find C go as a function of C ow : The horizontal asymptote for the above equation is which corresponds to the situation where all the water has been displaced from the corner.
The simulation on the coarse, uniform grid uses grid-cell spacing x = 1, reinitialization frequency once every 5 MLS iterations, and c = 0.001 as MLS convergence criterion in equation (18).Now, halving the grid spacing will reduce the CFL-determined time-step four times.So, to make simulations comparable between meshes with different grid spacing, we use a reinitialization frequency of once every 20 iterations for AMR-1 and once every 80 iterations for AMR-2.To maintain comparable convergence criteria, we also set c = 0.0005 for AMR-1 and c = 0.00025 for AMR-2 in equation (18).Simulations on the corresponding uniform, fine grids with x = 0.5 and x = 0.25 follow the same procedure.In the AMR simulations we tag cells for refinement around interfaces in the pore space according to the criterion |φ α | < x and ψ > 0 for α = g, o, w.We perform regridding every 1000 iterations for AMR-1 and every 2000 iterations for AMR-2.
Figure 4 shows three-phase fluid configurations for P gw = 10 and 20 kPa from the simulation on the finest uniform grid ( x = 0.25) and the AMR-1 and AMR-2 simulations.As the AMR simulations form grids with higher resolution only around the fluid/fluid interfaces, the pore geometry is resolved with coarse grid spacing ( x = 1) in areas away from the interfaces, leading to rounded pore corners in case of AMR-1 and AMR-2.Such differences in uniform (or initial) grid spacings also lead to minor differences in the input interface location on the computational grids, which contribute to slightly different original volume of the oil layers in the corners.This difference is larger for the 30 • corner in comparison to the 60 • corner.The supplementary information provides zoomed-in views of the configurations in the different pore corners for P gw = 10, 14, and 20 kPa.
A comparison of the results from the AMR-1 simulation with a simulation on the corresponding uniform, fine grid ( x = 0.5) shows that the locations of the zero contours of φ α from the two simulations are nearly overlaying each other, see Fig. 5 (a) and (b).The deviation between the two simulations is largest for the oil/water interface in the 30 • -corner where the zero contours are separated by at most one coarse grid spacing ( x = 1).This separation can be attributed to the slightly different initial oil-layer volumes, which leads to a more pronounced separation as the oil is squeezed deeper into the corner.The overlapping C go (C ow )-curves shown in Fig. 5 (c) and (d) confirm that the AMR simulation is able to mimic the results expected from the corresponding uniform fine grid.
Figure 6 presents a similar comparison between the AMR-2 simulation and the simulation on the corresponding uniform, fine grid with grid spacing x = 0.25.As for the previous case, the zero contours of φ α are nearly overlapping, Figure 7 presents C go (C ow )-curves from the AMR-1 and AMR-2 simulations along with the results from the simulation on the coarse uniform grid ( x = 1).Clearly, the coarse-grid simulation reasonably resolves the interface curvatures in the 30 • -corner, while AMR improves accuracy of the solution further (see Fig. 7 (a)).The resolution of oil-water interface is more accurate in the 30 • -corner compared to the 60 • -corner.So, the most significant improvements with AMR occurs in the 60 • -corner, where the results on the coarse, uniform grid noticeably differs from the analytic solution (see Fig. 7 (b)).These results support the application of AMR in complex 3D geometry to improve the accuracy of results.Figure 8 explores convergence behavior of the oillayer volumes based on these three simulations.We present results from the equilibrium states obtained at every 2 kPa increase of P gw and estimate the relative error E[V ] of the oil volumes calculated during the MLS iterations with respect to the initial (target) volume in each corner.In the case of the 90 • corner, the relative error for different P gw almost coincide because the absence of water renders the gas-oil interface static.All results show that the application of AMR in the MLS-LVC model yields quadratic volume convergence with respect to grid spacing.The data is from uniform grid simulation ( x = 1), AMR-1, and AMR-2 simulation.The x-axis represents the grid spacing on the finest grid in the simulation Finally, Fig. 9 depicts various aspects of the computational performance of the MLS-LVC model with AMR.The total number of grid cells decreases during the AMR simulations because fewer grid cells are tagged for refinement when the interfaces approach the corners, see Fig. 9 (a) and (b).We also note that even the AMR-2 simulation forms a lower number of grid cells than the uniform-grid simulation with grid spacing x = 0.5.Similarly, Fig. 9 (c) shows that the numbers of patches decrease during the AMR simulations, even though they remain higher than the constant patch numbers from the simulations with uniform grid.
Figure 9 (d) compares CPU time spent per 1000 iterations in the simulations.The coarse-grid simulation ( x = 1) is fastest, while the AMR simulations are faster than the corresponding uniform fine-grid simulations.The number of iterations required to reach the final stable state is the same for the corresponding simulations with AMR and uniform fine grid.All the simulations become progressively faster as the oil layers move deeper into the corners.This is mainly because the periodic reinitializations of φ α require less iterations in the late stage of these simulations.In gen-eral, a dynamic CPU time will also arise when the number of patches spanned by conserved domains changes significantly during a simulation, leading to varying communication loads between processors in the LVC method.Specific to the AMR simulations on the triangular pore is that both the numbers of grid cells and patches decrease, although a generally higher number of patches in the cases with AMR can maintain a significant processor communication in the LVC method and lead to more modest declines of CPU time.Still, Fig. 9 (d) shows that, on an average AMR-1 is twice as fast, and AMR-2 more than four times as fast, as the computation on a corresponding uniform fine grid.Combining this efficiency gain with the nearly overlapping results from Fig. 5 (c) and (d) and Fig. 6 (c) and (d), we conclude that AMR is computationally important for higher accuracy in large geometries.

Validation against three-phase fluid configurations in a 3D cylindrical pore
We consider three-phase fluid configurations in a 3D cylindrical pore where water is the wetting phase, oil is the intermediate-wetting phase, and gas is the non-wetting phase.
Water is present at the two ends, while two oil droplets of different sizes surround a single gas bubble (see Fig. 10 (a)).
We conserve the volume of the gas bubble and the two oil droplets, while we treat the water as a continuous phase with a prescribed pressure.We provide the initial configuration as cylindrical regions (fluid/fluid interfaces are parallel to x-y plane) and allow the gas-oil and oil-water curvatures to develop during simulation based on the fluid phase properties (see Fig. 10 (a)).For a given water pressure P w , the oil pressure P o , and gas pressure P g can be calculated analytically from Young-Laplace equation as: r , where r is the pore radius.Note that two oil droplets with different (but sufficiently large) volumes will have identical pressures in such straight tubes.
For the validation, we study the fluid configurations in three pores with length 90 μm, and pore radii 10 μm, 12 μm, and 15 μm.The spatial domain size is 40L * × 40L * × 90L * .The characteristic parameters for the MLS-LVC method are L * = 10 −6 m, γ * = 10 −4 Nm -1 , and P * = 100 Pa.Hence, the coarse grid spacing is x = 1.We model a slightly non-spreading fluid system with interfacial tensions σ ow = 20 × 10 −3 Nm -1 , σ go = 11 × 10 −3 Nm -1 and σ gw = 30 × 10 −3 Nm -1 .The contact angles are θ ow = 20 • , θ go = 24.24• and θ gw = 16.1 • .We prescribe the water drop pressures to 2 kPa.All other simulation parameters are the same as in the last section.During simulation, the pressures of the two oil droplets and the gas bubble will be calculated from equation (20).Figure 10 (b-d) shows the equilibrium fluid interfaces for pore radius 15 μm.The supplementary information contains images of the equilibrium fluid configurations for the other pore radii.
Table 1 lists the pressure and volume of the gas bubble and the two oil droplets from our simulations and their analytical values.Figure 11 compares analytical and simulated gas and oil pressures, while Fig. 12 shows the convergence of the gas and oil volumes with respect to the grid spacing.Figure 12 confirms a quadratic volume convergence for the Fig. 11 Comparison of simulation results vs. analytical solutions for oil droplet and gas bubble pressures in the three cylindrical pores of different sizes.In each simulation, the two oil droplet pressures are nearly identical and appear as one data point in the figure.The three leftmost data sets represent oil droplet pressures, and the three rightmost data sets represent gas bubble pressures different pore sizes.Table 2 lists the computational data for these simulations.It can be seen that the number of iterations required for convergence increases with pore size.The computational time is also affected by the changes in the number of grid cells and patches in the geometry (computational performance analysis is discussed in detail in the next section).The results show that increasing the refinement level improves the accuracy of the pressure and volume results.

Scaling tests
In this section, we will discuss the parallel performance and scalability of the MLS-LVC method with AMR.For uniform grids, Jettestuen et al. [28] demonstrated that the algorithm performance depends on the number of preserved domains.Here, we investigate the parallel performance with one and two levels of AMR, which are the most relevant cases for application of the method on segmented micro-CT images of rock samples.We carry out weak scaling and strong scaling tests for uniform coarse grid, AMR-1, and AMR-2, while conserving domains of one phase (LVC-1) and two phases (LVC-2).
In the scaling tests, we consider a pore space occupied by an equal number of disconnected oil and gas domains (see Fig. 13).Continuous water occupies the remaining part of the computational domain (solid is absent).Pressures of continuous phases are set equal, while pressures of disconnected domains are calculated from equation (20).The simulations consist of 24 level set iterations for each of the three fluid phases, including the projection step calculations as well as the LVC calculations for selected phase(s).Three reinitialization solves (each consisting of iterations) are carried out at the 10th and 20th level set iteration step.In AMR simulations we refine all grid cells that satisfy |φ α | < x, α = g, o, w.The simulation time in the scaling tests is the duration Fig. 12 Volume convergence of oil droplets and gas bubble with respect to the grid spacing in the cylindrical pore: (a) r = 10μm, (b) r = 12μm, and (c) r = 15μm.E[V ] is the relative error, and the x-axis represents the grid spacing on the finest grid in the simulation.Gas in green-solid, smaller oil drop in red-solid, and larger oil drop in blue-dashed between the 4th iteration step and the 24th iteration step and thus we exclude gridding computation times.We refer to Gunney and Anderson [14] for an investigation of the parallel performance of AMR regridding within the SAMRAI framework.Even though the scaling performance of SAM-RAI is linear for processor clusters much larger than the 512 processors used in our test, the scaling test performance presented in this section should diminish with the frequency of regridding during the test.Overall, the scaling tests capture 60(= 20 × 3) level set iteration steps and 60(= 2 × 10 × 3) reinitialization steps.All the computations are performed on the supercomputer Fram provided by UNINETT Sigma2the National Infrastructure for High-Performance Computing and Data Storage in Norway.
In the case of strong scaling, we consider a computational domain with 320 × 320 × 320 coarse grid cells for the uniform-grid and AMR-1 simulations (see Fig. 13 (a)).As this domain size leads to memory overflow for AMR-2 with less than 512 processors, we use a computational domain with 320 × 320 × 160 coarse grid cells for the strong scaling tests with AMR-2 (see Fig. 13 (b)).In both cases the geometry contains a total of 216 (= 6 3 ) equally sized isolated fluid domains.Figure 14 shows the results from the strong scaling tests, and Tables 3 and 4 detail the speed-up and number of patches data.The tests show that AMR has a limited effect on the speed-up behaviour of the MLS-LVC method.In general, LVC-1 provides better scalability compared to LVC-2.Compared to the results on the uniform coarse grid, the scalability dips for AMR-1 but becomes equal or better for AMR-2.The number of patches per processor are nearly constant for the simulations on the uniform coarse grid, whereas in the AMR runs, the number of patches per processor decreases as the number of processors increases.The non-monotonic performance is due to a more significant decrease in the number of patches per processor for AMR-2 compared to AMR-1.Such a decrease reduces the load on each processor, and it also reduces the communication requirement and root processor computations involved in identifying unique isolated domains in the LVC algorithm.
In case of weak scaling, we consider a cubic domain with 320 × 320 × 320 coarse grid cells for 512 processors and 64(= 4 3 ) isolated fluid domains (see Fig. 13 (c)).Then we proportionately reduce the number of coarse grid cells and the number of processors by a factor of eight until we obtain 40× 40 × 40 coarse grid cells for one processor, while the number of isolated fluid domains (64) remains constant.Tables 5 to 7 lists the results of the weak scaling tests.For the smallest and least resolved geometry (using one processor), the first level of refinement in the AMR simulations covers the complete geometry instead of being limited to the region around the interfaces.This leads to a longer computational time than for the AMR simulations on the larger geometries with 8, 64, and 512 processors.Another trend is that the computational time for the AMR runs is least for 64 processors, see Table 5.This trend results from two parameters: the number of patches per processor and the total number of patches.Like the strong scaling test, the reduction in the number of patches per processor can increase computational speed.However, an increase in the total number of patches increases communication time and computational time on the root processor.We illustrate the impact of the total number of patches by comparing computation times for AMR-2 runs with 64 and 512 processors (Table 5) against corresponding patch numbers in Table 6: Even though the number of patches per processor is significantly lower for the 512-processor runs than for the 64-processor runs, the 512-processor runs are slower as they form a nearly five times higher number of patches in total.These findings demonstrate that weak scaling results are sensitive to the total number of patches in the system, as reported by Gunney and Anderson [14].The SAMRAI framework has been designed such that, for given grid refinement criteria, the scalability of the regridding process is of a linear order, but there is no control on the total number of patches that will be formed in the system.
An outlier in the weak scaling results is the 512-processor runs for AMR-1 which is slower than the 8-processor runs.In this case, the number of patches per processor has increased in comparison to the 64-processor case.The number of patches for 512 processors is approximately 24 times higher than for 8 processors.Based on the reasons mentioned above, there is a compounded decrease in performance which makes the overall computation slower.Table 7 highlights the reduction in the number of grid cells in AMR computations compared to a corresponding uniform fine grid.These values confirm the importance of using AMR to reduce the size of a large computational grid and using a large number of processors.Additional results from the scaling tests are provided in the supplementary information.The SAMRAI framework provides flexibility in patch generation by taking user inputs for the largest and smallest patch sizes during simulation.In our model, we use the computational domain size as the largest patch size and twice our numerical stencil size on each axis as the smallest patch size.This choice allows us to prevent a large number of small patches that can hinder computational speed.An important factor impacting the computational speed is the number of isolated ganglia in the system, which is apparent by the longer computation time of LVC-2 simulations compared to LVC-1 simulations.

3D LS/MLS-LVC simulations with AMR on sandstone pore geometries
We use 3D pore geometries of sandstone to explore the impact of AMR in simulations of capillary-controlled twophase and three-phase displacements where disconnected fluid ganglia form.For this purpose, we have extracted subvolumes from a data set of segmented microtomography images of Castlegate sandstone with voxel length 5.6 μm [58].We use LS-LVC method for two-phase simulations and MLS-LVC method for three-phase simulations.In both cases we compare the results produced with one level of refinement (AMR-1) against the results from the simulation on the coarse, uniform grid (represented by voxels).We set the MLS characteristic parameters for the simulations to L * = 5.6 × 10 −6 m, γ * = 10 −4 Nm -1 , and P * = 17.9 Pa.The number of conserved domains doubles in LVC-2 simulations compared to LVC-1 simulations.For example, in the weak scaling test, there are 32 and 64 locally conserved domains in LVC-1 and LVC-2 simulations, respectively.In contrast, the number of patches per processor typically remains equal since we refine grid cells around all fluid/fluid interfaces in the pore space (see Tables 4 and 6).The results from the strong and weak scaling tests with 512 processors present four different numbers of conserved domains on the 320×320×320 geometry.The results show that the factor by which the computational time increases with a growing number of locally conserved domains is lower than the growth in the number of such domains (see supplementary information for details).

Two-phase imbibition
We explore two-phase oil/water displacements during imbibition on a Castlegate rock geometry with volume 200L * × 200L * × 200L * and porosity 22.7%.The rock is considered to be water-wet with contact angle θ ow = 20 • , and interfacial tension is σ ow = 20 × 10 −3 Nm -1 .For this case, Helland et al. [17] performed LS simulations of primary drainage and imbibition without LVC, by increasing (drainage) and then decreasing (imbibition) the interface cur-vature C ow = σ ow /P ow stepwise after each equilibrium state achieved with equation (7).Here, we simulate imbibition with LVC applied to oil, for both AMR-1 and coarse uniform grid, starting from a reversal point on the primary drainage P ow (S w )-curve where S w = 0.13 and P ow = 2.68 kPa.We treat the top face of the rock as inlet boundary where we place a water-wet porous plate with thickness of two coarse grid cells, which allows water to imbibe into the sample.Oil withdrawal occurs only through bottom face (outlet) while the side walls are sealed.The level set functions use reflective conditions at all boundaries.We terminate the imbibition simulations when all the oil is capillary trapped as disconnected ganglia.The coarse, uniform simulation uses grid spacing x = 1, while periodic reinitialization of φ occurs every fifth LS iteration.The tolerance value in equation ( 18) is c = 1.5 × 10 −3 .Based on the discussions in section 4.1, we let the AMR simulation perform reinitializations every 20th LS iteration and set c = 7.5 × 10 −4 .We tag cells for refinement that satisfy |φ| < x, ψ > 0 and regrid every 2000 iterations.
Figure 15 shows the location of oil phase and finer grid patches at three different capillary pressures during the AMR Fig. 16 Residual oil configurations after imbibition in Castlegate sandstone from AMR-1 and uniform, coarse-grid simulations.Images (a) and (b) show oil-phase locations from two diagonally opposite views for AMR-1 (orange), coarse uniform grid (cyan), and the common locations from both simulations (red) Fig. 17 Capillary pressure between continuous oil and water as a function of water saturation for primary drainage (coarse, uniform grid) and imbibition (coarse, uniform grid and AMR-1).Capillary pressures of disconnected oil ganglia that form during imbibition for AMR-1 (orange circles) and coarse, uniform grid (purple crosses) are also shown simulation.A video of the AMR simulation is provided as supplementary information.During the initial stages of imbibition, the oil phase is continuous while water layers in crevices along the pore walls grow in thickness [17].Even- tually, water snaps off oil and invades narrow pore throats while disconnected oil ganglia form and relax to stable states in nearby pore openings.Figure 16 shows that the location of oil ganglia changes slightly due to refinement.However, Fig. 19 Three-phase fluid configurations of water (blue), oil (red) and gas (green) from AMR simulation of gas invasion in Castlegate sandstone for (a) P gw = 1.9 kPa (initial state), (b) P gw = 2.77 kPa, and (c) P gw = 3.44 kPa (final state).Patches with fine grid are shown in black the residual oil saturation after imbibition is S or = 0.27 for both the uniform coarse grid and AMR-1 runs.
Figure 17 reveals that the capillary pressure curves from the two simulations are almost identical and that most oil ganglia form when S w > 0.4.An important point of note is that interface relaxation creates stable states in which the capillary pressures of oil ganglia differ from the prescribed capillary pressure between the continuous phases at the incident of ganglia disconnection.The oil ganglia create a wide range of capillary pressures, in which the majority of the data are located between the continuous P ow (S w )-curves for drainage and imbibition.Such behaviour are also observed in pore-scale experiments based on analysis of two-phase micro-CT images [12,20,34].The AMR simulation differentiates the ganglia capillary pressures more distinctly as it creates a wider spread among the majority of the data than the coarse simulation.
Figure 18 (a) shows that the number of patches remains nearly constant while the refined area increases during the initial stages of the imbibition.The increase in refined area is Fig. 20 Simulation results from gas invasion on Castlegate sandstone using AMR and coarse, uniform grid: (a) Saturation paths, and (b) gaswater capillary pressure curves due to the swelling of water layers which increases interfacial area between oil and water [11,17].At later stages, both the number of grid cells and the number of patches in the computational geometry decrease as oil is displaced.This speeds up the simulation and results in decreased slope of the total CPU usage curve, as shown in Fig. 18 (b).The total CPU time is the product of simulation time and number of computing cores.
In the present study, we repeat the simulation of gas invasion with one level of grid refinement (AMR-1).The inlet is located at the bottom face of the cubic domain where we add a two-voxel thick layer of pore space filled with invading phase.With respect to LVC, the side walls are sealed and top face is outlet.The level set functions use linear extrapolation at inlet and reflective conditions at the other faces as boundary conditions.We initiate gas invasion from a high water saturation (S w = 0.67) before oil breakthrough occurs in primary drainage.Thus, all the oil is disconnected and subjected to LVC in the beginning of the three-phase displacement.We simulate gas invasion by increasing gas pressure stepwise while the phase pressures of continuous oil and water are held constant.The simulations are terminated when the oil saturation no longer decreases.The coarse-grid simulation performed reinitialization of level set functions φ α every 10th MLS iteration and used c = 1.5×10 −3 as tolerance value.As we are mainly interested in the behaviour of the disconnected phase in the AMR simulation, we only refine cells around oil boundaries according to |φ o | < x and ψ > 0. Gas/water interfaces will then be located on the coarse grid, so we maintain a reinitialization frequency of once every 10 iterations while we use a reduced tolerance value c = 7.5 × 10 −4 as discussed previously.
Figure 19 shows three-phase fluid configurations and patches with finer grid for three different gas-water capillary pressures during the AMR simulation.A video is provided as supplementary information.In the initial stages, gas invasion occurs by double displacements in which gas displaces several oil slugs that displace water.After oil reaches the top boundary the dominating mechanisms are direct gas-oil  [28] displacements (that reduces oil saturation) and double displacements.During this process, oil ganglia detach from both continuous oil and larger oil slugs of which some of the ganglia continue their displacement while others get trapped.The trapped oil phase occupies mostly intermediate-sized pores and pore corners as intermediate-wet layers.The last stage of the simulation is dominated by direct gas-water displacements.
Figure 20 compares the three-phase saturation paths and gas-water capillary pressure curves from the AMR-1 and coarse-grid simulations.Typically, capillary pressure curves for drainage on small pore-scale samples display a staircase trend (for pressure-controlled invasion) that consists of intervals with irreversible saturation jumps between equilibrium states (e.g., Haines jumps) connected by intervals with smooth reversible saturation changes [7,17,19].Although the level of the two capillary pressure curves are similar, the AMR simulation shows a more pronounced staircase trend than the coarse simulation.The saturation paths from the two simulations differ more.The deviation at intermediate gas saturations occurs because the AMR simulation creates a different route for the gas and oil movement to the top face, leading to gas breakthrough at a higher gas saturation.Figure 21 shows that this leads to significantly different residual oil configurations in which the AMR simulation traps more oil in the lower half of the geometry than the coarse simulation.The residual oil saturations from the simulations are S or = 0.156 (coarse grid) and S or = 0.182 (AMR).Further, the residual oil configuration from the AMR simulation is less fragmented as it contains fewer but more elongated oil ganglia with shapes that are more reminiscent of intermediate-wet oil layers from three-phase pore-scale experiments in porous rocks [e.g., 53]. Figure 21 shows that there is also some differences in the final gas configurations from the two simulations even though the final gas saturations are about the same.
Figure 22 illustrates the performance of the MLS-LVC model during the AMR simulation.The numbers of patches and fine grid cells in the geometry increases as gas displaces and traps oil.After gas breakthrough, these numbers reach a maximum and experience a small decline when direct gaswater displacements start dominating.The slope of the total CPU usage curve responds to the change in the numbers of patches and grid cells during the simulation.

Conclusion
This work has added structured, adaptive mesh refinement capabilities to locally, conservative level set methods for simulating two-phase and three-phase capillary-controlled displacement with disconnected ganglia at the pore scale.The AMR implementation of these LS/MLS-LVC methods was carried out within the patch-based software framework SAMRAI.We created new integer operators that are essential in preserving the local volume of isolated domains.The implementation also improved memory requirements and computations by preventing multiple identifications of any particular set of neighbouring domains.We validated the MLS-LVC method with AMR against analytical three-phase configurations in idealized pores, where we treated the oil layers as distinct domains with their volumes preserved.The simulation results provide pressure and volume convergence while tracking more than one isolated domain.Specifically, we obtain second-order volume convergence with respect to grid spacing.The computational performance of the model during simulation underscores the strength of AMR compared to a corresponding uniform fine grid.Further, the results from simulations with AMR and a corresponding fine uniform grid coincide, which shows that AMR simulations provide similar accuracy with faster computations.The AMR implementation shows good performance in both strong and weak scaling tests.The scaling tests also point out the impor-123 tance of the number of patches per processor as well as the total number of patches in gaining computational efficiency.
Finally,t we investigated the impact of AMR on LS/MLS-LVC simulations of two-phase imbibition and three-phase gas invasion with formation of disconnected oil ganglia on 3D pore structures of water-wet Castlegate sandstone.Compared with coarse-grid simulations, the most significant impact of AMR was obtained for three-phase displacements with respect to the behaviour of the conserved, intermediate-wet phase.Differences in the shapes and spatial distribution of oil ganglia lead to different residual oil saturations and saturation paths.These results suggest that AMR is an important technology for obtaining desired accuracy and efficiency in pore-scale simulations of threephase displacements.This will in turn, improve calculations of constitutive relationships, such as three-phase capillary pressure-saturation curves, which are required for large-scale multiphase flow simulations in porous media.Although not intuitive, Table 8 shows that the parallel identification does not provide any gains.Instead, it is generally slightly slower than serial identification.The drop in speed can be attributed to three factors.First parallel identification requires an extra iteration of reading and writing data on a patch along with a processor level run for the unique labelling of domains.These additional executions increase the time spent on each processor before unique global labelling.Second, over the years, there has been a continuous reduction in MPI data communication time between processors, and nowadays, it is of linear order [14].This reduction diminishes the possible gain from a reduction in computational time from parallel processing.Finally, in our model we utilise the patch generation and load balancing The table lists the time for computation, t, the speed-up with parallel identification, S, the reduction in the number of elements compared to a uniform grid, E, and the number of patches, Pat, for a particular number of isolated domains, N d capabilities of the SAMRAI framework.SAMRAI provides adequate load balancing by striving to put nearby patches into the same processor.However, using the framework reduces our control in forcing adequate overlap onto the patches in a processor to ensure that parallel identification improves computation time.

Fig. 1
Fig.1AMR grid structure in SAMRAI framework with three patch levels.The regions tagged for refinement are shown in green.The illustration shows a case where each tagged grid cell is refined into four equally-sized grid cells.The calculations at any spatial location are done on the highest available patch level and then transferred to lower levels through ghost cells.For details, see Anderson et al.[3] (a) Isolated domain identification (Algorithm 1, Section 3.2): Construct domains { a : φ α ≤ 0 and a − b ∞ > 1, ∀ a = b, a and b ∈ I α (t 0 )}.Here, a − b ∞ represents the Chebyshev distance between the boundaries of domains a and b .Figure 2 (top row) depicts isolated domains identified under different levels of mesh refinement.(b) Generate a set of extended domains { ext a : a ∈ I α (t 0 )} (Algorithm 2, Section 3.2). Figure 2 (middle row) shows the result after extension with different levels of refinement.(c) Modify extended domains (Algorithm 3, Section 3.2): Modify the set of extended regions { ext a : a ∈ I α (t 0 )} by extending the domain labels across patches.The modification ensures that the numerical stencils for grid cells near patch borders and fluid/fluid interfaces use the correct pressures in the numerical scheme.Figure 2 (bottom row) shows the result after modification of extension maps with different levels of refinement.(d) Set t n = t 0 :

Fig. 3
Fig. 3 Integer operators where x and y are integer values and INT_MIN is the minimum integer limit in C++.(a) Refining operator, (b) coarsening operator 1 for equal values, and (c) coarsening operator 1 for different values

Fig. 9 Fig. 10
Fig. 9 Performance of simulations with AMR and uniform grids as function of the number of MLS iteration steps.(a) Ratio between the numbers of grid cells for AMR and corresponding uniform, fine grid,

Fig. 15
Fig. 15 Location of oil phase (red) during AMR simulation of imbibition in Castlegate sandstone for (a) S w = 0.13 and P ow = 2.68 kPa (initial state), (b) S w = 0.39 and P ow = 1.07 kPa, and (c) S w = 0.73 and P ow = −0.54kPa (final state).Patches with finer grid are shown in light blue

Fig. 18
Fig. 18 Model performance during the AMR-1 simulation of two-phase imbibition on Castlegate sandstone.(a) Evolution of the number of patches and the ratio of the number of grid cells between AMR-1 and a corresponding uniform fine grid.(b) Total CPU usage curve

Fig. 21
Fig. 21 Oil (top) and gas (bottom) configurations at the final state (P gw = 3.44 kPa) of gas invasion in Castlegate sandstone from three-phase simulations with (a) AMR-1 (this work), and (b) uniform, coarse grid[28]

Fig. 22
Fig. 22 Model performance during three-phase AMR-1 simulation on Castlegate sandstone.(a) Evolution of the number of patches and the number of grid cells (normalised by the number of grid cells in a corresponding uniform fine grid).(b) Total CPU usage curve

Fig. 23
Fig. 23 Geometry used for strong scaling tests with isolated domains of one phase (in blue) surrounded by a continuous phase that occupies the remaining part of the computational domain.(a) Geometry used for 16, 32, and 64 processors contains 50 × 50 × 50 coarse grid cells.(b) Geometry used for 128, 256 and 512 processors contains 160×160×160 coarse grid cells 1 , { a (t 1 ) : a ∈ I α (t 1 )} and { ext a (t 1 ) : a ∈ I α (t 1 )}, and similarly, at time t 2 = t 1 + t these sets are { b (t 2 ) : b ∈ I α (t 2 )} and { ext b (t 2 ) : b ∈ I α (t 2 )}.The intersection between the domain a (t 1 ) and the extended domain ext b (t 2 ) is ) 1: Initialise patch variable M to zero.2: Create an empty list L l on each processor where l is the rank of processor in MPI_COMM_WORLD.Also, create an empty list L agg on root processor.3: Update P to the rank of processor.Then, update P and M on the refined regions of a coarse patch to INT_MIN 4: Identify domains on a patch where φ α ≤ 0 using grassfire algorithm.sequentially to M global 26: Send relevant M global to each processor 27: Update M on the patches to M global Local ≥ 0) and (P Neighbour ≥ 0) and (M Local > 0) and (M Neighbour > 0) then 10: if (P Local < P Neighbour ) or (M Local = INT_MAX) then 11: Save neighbours tuple N in list L l 12: else if (P Local = P Neighbour ) and (M Local < M Neighbour ) agg ⊃ L l , ∀ l ∈ MPI_COMM_WORLD 20: On root processor, ∀ N ∈ L agg , use P as an indicator to reidentify patch-level unique M to globally sequential values M * .21: while on root processor and there exists N ∈ L agg : M * * value changed in the last step, update all the corresponding M * in L agg 24: end while 25: Reidentify all unique M *

end for Algorithm 4
Restructure variables during regridding Require: Patch variables: Phase level set value (φ α ), Processor rank (P), domain label (M), extended domain label (E), distance (D), ratio of patch grid spacing w.r.t.finest patch grid spacing (R), modification distance threshold (d c , in city-blocks), and number of grid cells in z-direction (n z ), y-direction (n y ) and x-direction (n x ).Array variables: original volume (V boundar y 1 11: if (d < d c and D > D ghostcell + d) or E i, j,k = 0 then 12: E i, j,k = E ghostcell 13: end if 14:

Table 2
Simulation performance in cylindrical pore for different pore sizes.Time in CPU hours

Table 8
Performance comparison between serial and parallel identification of isolated domains for a given refinement ratio, R, and number of processors, P.