A locally conservative multiphase level set method for capillary-controlled displacements in porous media

Article history: Received 28 January 2020 Received in revised form 27 October 2020 Accepted 28 October 2020 Available online 31 October 2020

We present a multiphase level set method with local volume conservation for capillarycontrolled displacement in porous structures. Standard numerical formulations of the level set method for capillary-controlled (or, curvature-driven) motions assume phase pressures and interface properties are spatially uniform and disregard the fact that separate phase ganglia typically have distinct pressures. This is a major problem for the suitability of such methods to simulate capillary trapping in porous rocks as it will lead to severe mass loss. The method presented here preserves volumes of individual phase ganglia, while it predicts capillary pressures between ganglia and surrounding phases. A conservative volume redistribution algorithm handles ganglia breakup and coalescence. The method distinguishes between three-phase systems, where separate level set functions describe the different phases, and two-phase systems, where one level set function represents interfaces. We present sequential and parallel algorithms for the new method and emphasize important aspects specific to the patch-based parallel implementation. We validate the method numerically by applying local volume conservation to simulations of two and three phase systems in both two and three spatial dimensions. The model is tested for both saturation and pressure controlled systems and handles both merging and splitting of phase ganglia.

Introduction
Understanding the mechanisms leading to fluid phase trapping and mobilization in three-phase flow in porous media are essential to strategies for improved oil recovery and carbon dioxide (CO 2 ) storage in mature oil reservoirs. This requires knowledge of how fluid ganglia behave, with coalescence and splitting, in the presence of other phases in the reservoir, as well as their impact on capillary pressure and relative permeability relations used in reservoir simulation. In most of the reservoir, away from wells and fractures, capillary forces typically dominate over viscous and gravity forces during multiphase displacements at the pore scale [1]. Therefore, ganglia of the non-wetting phase (e.g., gas or oil), surrounded by a continuous wetting phase (e.g., water), likely get trapped by capillary forces [2]. This is a desired mechanism for safe CO 2 storage [3]. However, in a three-phase system, the action of capillary pressure can mobilize ganglia by double displacement, in which a continuous phase (e.g., gas) displaces an isolated cluster of the second phase (e.g., oil) that in turn displaces the continuous third phase (e.g., water) [4,5]. A complexity with such three-phase displacements is their different behaviours in the presence of triple lines, where all three phases meet, and in the presence of spreading oil layers separating the other two phases. Pore-and core-scale experiments have shown that spreading oil layers, which make the oil phase more continuous, and double displacements, can mobilize oil and lead to low residual oil saturation after gas invasion [6][7][8]. An extension of the double-displacement mechanism refers to the situation where the invading phase displaces a chain of several neighbouring fluid clusters that displaces another continuous phase. This can also happen in two-phase systems under weakly-wet states when the fluid phase present on opposite sides of a ganglion are disconnected from each other and associated with different pressures. Multiple displacements play an important role as mechanisms for oil recovery over multiple water-alternate-gas injection cycles [9][10][11]8,12].
X-ray microtomography imaging of fluid distributions during multiphase flow experiments in porous media has allowed the analysis of shape and topology of fluid ganglia [7,12], capillary pressure states of ganglia compared with continuous phases [13,14], wetting state and displacement mechanisms [15,16,12]. The use of digital rock physics, specifically direct numerical simulation on segmented rock images, is potentially a tool for investigating three-phase processes. Fundamental to simulating multiphase flow with isolated ganglia on realistic pore structures is the numerical method's ability to conserve mass (or volume, for incompressible fluid flow) in a correct manner.
There are many methods for conservation of volume in multiphase flow. Lattice Boltzmann [17] and phase field [18,19] methods add interaction between phases and use the emerging phase separation as the basis for multiphase solvers, where the phase interaction potential set surface tensions and contact angles. In these schemes fluid masses are explicitly conserved but undesired numerical effects such as mutual pollution of non miscible phases, large diffuse interface regions, and spurious currents, can occur.
For capillary-dominated multiphase flow regimes, a reasonable assumption is that capillary forces alone control the interface displacement. Examples of such quasi-static, capillary-controlled methods include pore-network modelling that typically uses an idealization of the pore geometry [e.g., 20] and pore-morphology methods [21]. Intrinsically, these methods require additional rules or techniques for volume conservation.
Implicit interface-tracking techniques are direct numerical simulation methods suitable to capillary-controlled (or, curvature-driven) displacements. Within this category, the two most widely used approaches are the level set (LS) method [22][23][24][25] and the volume of fluid (VOF) method [26]. VOF methods have the advantage that they conserve mass but calculation of interfacial properties like normal vectors and curvatures are cumbersome. For the level set method it is the other way around. Calculating interfacial properties is straightforward, but the method has known issues with mass conservation. However, there are methods for improving mass conservation in LS methods. For instance, Sussman and Puckett [27] developed a hybrid method using both VOF and LS approaches, Enright et al. [28] used a hybrid particle-LS method to preserve volume, and Luo et al. devised a method based on spectrally refined interface approaches [29]. Saye and Sethian [30] preserved volumes of each phase within their LS-based Voronoi implicit interface method, by modifying the phase pressures to control interface motion in normal direction.
In this work, we develop a multiphase level set method with local volume conservation (MLS-LVC) to study capillarycontrolled three-phase displacements at the pore scale. We have previously implemented global volume conservation in our LS [31] and MLS [32] methods using the approach of Saye and Sethian [30], which allows simulations with prediction of capillary pressure for specified phase saturations [33,32,34]. Here, we extend this approach further to conserve the individual volumes of all isolated ganglia of one or more phases, while treating ganglia coalescence and splitting through a new conservative volume redistribution algorithm. The resulting MLS-LVC method predicts surface area, volume, and capillary pressure for individual fluid ganglia at equilibrium states, thus providing complementary insights to pore-scale experiments addressing ganglion behaviour [e.g., 14,7].
The paper is organized as follows: First, Section 2 presents the main features of the MLS method. Then, we describe the LVC method in the MLS setting (Section 3) and explain the modifications needed to implement it for two-phase systems using a standard LS approach (Section 4). Section 5 presents important aspects of our parallel implementation using LVC. Section 6 presents three-phase MLS-LVC simulations of double and multiple displacements through a 2-D sinusoidal pore channel where several ganglia coalesce and split. We investigate both the effects of applying LVC to one and two phases as well as the behaviour of volume conservation and capillary pressure convergence of disconnected phase ganglia as grid resolution increases. Thereafter, we apply the LVC algorithm to the LS and MLS methods to carry out respective two-and three-phase simulations of primary drainage and subsequent gas invasions in a synthetic 2-D porous medium and on a small subvolume of 3-D sandstone. Section 7 summarizes our work with main conclusions. Future work will target simulation of three-phase residual saturations after water-alternate-gas invasion cycles for comparison with experiments and empirical correlations.

Multiphase level set (MLS) method
The level set method [22][23][24] is a numerical scheme for propagation of curves and surfaces in accordance with a given velocity function, F . This is accomplished by defining a curve or surface implicitly as the zero contour of a scalar function φ( x), where φ( x) is the position vector, that evolves according to the equation: (1) Fig. 1. MLS representation of fluid/fluid/solid contact lines in a pore throat containing fluid phase α and β before (left) and after (right) enforcing the projection step of Losasso et al. [35]. (For interpretation of the colours in the figures, the reader is referred to the web version of this article.) It is customary to split the velocity term into normal and advective components, F N and F A , respectively, which leads to the following equation: A strength of the level set method is that normal vector n of an interface and its curvature κ is calculated from derivatives of the level set function so that n = ∇φ/|∇φ| and κ = ∇ · n.
For systems containing more than two fluid phases, we will use a multiphase level set (MLS) approach [34,32] that assigns one level set function and one evolution equation per fluid phase α. This yields We identify domains with φ α < 0 as part of fluid phase α. The static solid phase is also defined by a scalar field ψ , where ψ < 0 is identified as part of the solid and ψ > 0 is pore space. One challenge with using multiple level set functions is the presence of overlap regions, where two or more phases are present, and vacuum regions, where no phases are present.
To avoid these issues we use a projection step method described by Losasso et al. [35], which amounts to identifying the two smallest level set functions at each point of the domain and then subtract the average of these two from all level set functions. We will investigate quasi-static, capillary-controlled multiphase systems [36,34,32], with well defined contact angles between fluid/fluid interfaces and pore walls [31,34]. The bulk part of the system is governed by the Young-Laplace equation for the fluid/fluid interfaces, where p α and p β are phase pressures, P αβ is the capillary pressure, σ αβ is the interfacial tension, and C αβ is the interface curvature, for interfaces between phases α and β.
The contact angle θ αβ , measured through phase β, is the angle between the αβ-interface and the solid pore wall. It is given by the fluid/fluid interfacial tension, σ αβ , and the solid/fluid interfacial tensions σ sα and σ sβ through Young's equation: The MLS method assigns a surface tension contribution γ α and a solid/fluid intersection angle β α to each fluid phase [34,32]. The γ α 's are related to the fluid-fluid interfacial tensions through the identity [36] σ αβ = γ α + γ β . (6) We define the solid-fluid interfacial tensions as σ sα = γ s − γ α cos β α , where γ s is chosen so that all σ sα > 0. Further, β α is a specified angle between the pore wall and the boundary of phase α, defined by cos β α = n α · n s , where n α and n s are the normal vectors derived from φ α and ψ , respectively, see Fig. 1. Consistent with Eq. (5), the contact angle θ αβ is given by cos θ αβ = n β · n s = γ β cos β β − γ α cos β α σ αβ . (7) Note that this expression is independent of γ s . For given sets of σ αβ and θ αβ we use γ -and β-parameters that satisfy Eqs. (6) and (7) as input to the MLS simulations. In two-phase simulations, using two level set functions, a possible choice for the intersection angles is β β = θ αβ and β α = 180 • − θ αβ .
The algorithm used to find the local equilibrium of a multiphase system, where each phase is given a uniform pressure, is based on assigning an energy cost γ α for each phase boundary so that the total energy of the system becomes the sum of bulk and surface contributions for each phase separately. We have explored this approach previously [32,34]. Here, we also want to consider that the pressure in each isolated region of each phase can be different, so that we need to keep track of individual regions that can change shape and topology in time. However, we shall assume that all regions of the same phase still have the same surface tension and contact angle. The total energy of the system [32] can now be written as where is the computational domain, H(· · · ) is the Heaviside step function and δ(· · · ) is the Dirac delta function. The velocity functions that are used to minimize this energy are [32] F α where κ α = ∇ · (∇φ α /|∇φ α |) is the curvature field calculated from phase α's level set function, and S(· · · ) is the sign function. We note that the pressure p α , and consequently κ α , is not only dependent on phase but also, possibly, on domain as will be explained in the next section. We use an explicit numerical scheme based on finite differences to solve Eqs. (3), (9) and (10), as described elsewhere [32,24]. With the projection method of Losasso et al. [35] enforced at the end of each iteration-time step, the equilibrium solutions satisfy [32] (p α − γ α κ α ) |∇φ α | = p β − γ β κ β |∇φ β | (11) in pore space, and γ α (∇φ α · ∇ψ − cos β α |∇φ α ||∇ψ|) = γ β ∇φ β · ∇ψ − cos β β |∇φ β ||∇ψ| (12) at the pore walls. Around αβ-interfaces, we have that φ α = −φ β , which implies ∇φ α = −∇φ β , |∇φ α | = |∇φ β | and κ α = −κ β . Thus, Eqs. (11) and (12) are consistent with Eqs. (4) and (7).

Local volume conservation (LVC) method
Our goal is to develop a simple and lightweight algorithm that, to a good approximation, conserves the volume of a phase and each isolated region of that phase. When evolving this system, regions can coalesce and split, so we need to devise an algorithm that can robustly redistribute volumes between regions that change both shape and topology. We introduce the following notation: For a given fluid phase α, let I α (t) ⊂ N be the enumeration of all isolated regions of that phase at iteration time t, similar to an indicator function. Then, a (t) ⊂ , where a ∈ I α (t), is an isolated region in space. To implement volume conservation of the isolated regions we use the method introduced by Saye and Sethian [30] where the pressure p a (t) of region a , a ∈ I α (t), at iteration-time t, is given as Here, V a (t) is the target volume of region a , V a (t) is the region volume at iteration-time t, A a (t) is the region surface area, and t is the iteration-time step. The pressure p a only equals the pressure in region a when the evolution of the LS fields has reached its equilibrium state. In this stationary state, V a will deviate from V a by a small amount governed by t. Hence, this method works by assigning an artificial numerical compressibility to every isolated region. We note that t is not a physical time so that its dimensions are chosen according to Eq. (13), also in agreement with the iterative time unit in Eq. (3).
We identify isolated regions with the grassfire algorithm described by Prodanović et al. [37]. To each region we assign a volume V (0) a so that the total target volume of phase α is given as Our method requires a strategy to assign volumes to regions at the next time based on the volume of regions at the previous time step. For this purpose, we define extended regions ext a (t) ⊃ a (t), a ∈ I α (t), that completely partition the computational domain, i.e., = ∪ a∈I α (t) ext a (t). We construct these extended regions ext a (t) by including all points nearest to a given region a (t) using a distance transform. In this work we have used the method described by Rosenfeld and Pfaltz [38]. Using this algorithm, we can calculate the integer distance d( x, a ) from a voxel (or, grid cell) at position x to the boundary of each region a , a ∈ I α . For each voxel x / ∈ a , a ∈ I α , we mark it as part of the extended region ext a , given The stippled lines show the position of interfaces in the previous/next time step. For the case with domain coalescence (t 2 > t 1 ), we obtain the following target volumes at t 2 : . For the case with domain splitting (t 1 > t 2 ), we obtain the following target volumes at t 1 (by swapping a and b, and t 1 and t 2 , in Eq. (16)): We define the intersection between region a at time t 1 and the extended region ext b at time t 2 > t 1 as Then, we use the relative volume of overlap between the t 1 and t 2 regions as a weight for the redistribution of volume. The target volume assigned to region b (t 2 ) is thus given as (16) where | · · · | denotes the volume of a region, measured in the number of voxels it contains. In our scheme we know that the ext regions cover the whole computational domain so that no volume is lost in the volume reassignment. Note that when coupled with MLS iteration scheme, restrictions on time step lead to relatively small changes in the fluid configurations within an iteration step. This ensures that target domain volumes are redistributed correctly with LVC. Fig. 2 illustrates the volume redistribution algorithm for simple examples with domain coalescence and splitting.
The pressure term given in Eq. (13) for region b (t 2 ), b ∈ I α (t 2 ), is obtained by using V (0) b (t 2 ) as the target volume and by using the following volume integrals to evaluate the volume, V b (t 2 ), and area, A b (t 2 ): H(ψ )δ(φ α )|∇φ α |dV . (18) Although the accuracy of these volume and area calculations can decrease when domains are sufficiently close to each other, typically on the length scale of the neighbourhood used by the numerical stencil, the sum of target volumes, V To allow for mass entering and/or leaving the system, we tag all regions of the same phase in contact with inlet/outlet boundaries of the computational domain. We assume all such regions are connected to the same fluid reservoir and belong to the continuous phase with the same pressure. This phase pressure is either assigned or calculated by controlling the phase saturation [32] (see also Section 6). Regions tagged as boundaries are excluded from the local conservation scheme as the pressure is given by the continuous phase. Once new isolated regions detach from the continuous phase, we calculate their target volumes with Eq. (17). Then, Eq. (13) implies that the phase pressures for such regions are zero in the first iteration after region detachment. However, in subsequent iterations these phase pressures eventually converge to correct values as the system relaxes toward equilibrium.
For simulations on realistic pore geometries it is inevitable that tiny phase domains containing only a few voxels arise. This can lead to high and unrealistic region phase pressures calculated from Eq. (13), which in turn can make the MLS iteration scheme unstable. In this work, we resolve this issue by setting the phase pressure of regions with target volumes smaller than a few grid-cell volumes equal to the pressure of the corresponding continuous phase. Specifically, we use 3( x) 2 (in 2-D) and 20( x) 3 (in 3-D) as volume limits. Although this implies that tiny regions can vanish during the MLS iteration steps, our simulations (see Section 6) indicate that the effect on the volume redistribution is negligible. We note that an alternative approach, that we have not investigated here, would be to exclude tiny regions completely from the LVC algorithm.

Formulation of the algorithm
We now give the main steps involved in the implementation of our LVC method in the MLS setting. AlgorithmSeq 1. Initialize the system: (i) Construct the set { a : a ∈ I α (t 0 )} (using grassfire algorithm) with unique tags a ∈ I α (t 0 ) assigned to each isolated fluid region of the conserved phase α, and (ii) generate the accompanying set of extended Loop over iteration time (from t n to t n+1 ): 2. Set pressures of conserved phase α equal to the continuous-phase pressure in regions connected to inlet or outlet boundaries. For the other isolated phase α regions, set pressures to p a , a ∈ I α (t n ), using Eq. (13), enforcing volume conservation. 3. Iterate the multiphase system from t n to t n+1 with Eqs. (3), (9) and (10), using the phase pressures from step 2, and then enforce the projection step [35].
4. Repartition the system at time step t n+1 , that is, construct the set Tag regions as boundaries if they belong to the phase connected to inlet or outlet boundaries.

Extend the partitioning
6. Redistribute volumes according to Eq. (16), using t n and t n+1 as time step t 1 and t 2 , respectively. Regions tagged as boundary belong to the continuous phase and are excluded from the volume redistribution. 7. Identify new regions that detached from the continuous phase and calculate their region target volumes based on Eq. (17).
After repeating this scheme a specified number of times, we reinitialize the φ α -functions to signed distance functions by solving [39] The iteration-loop ends when the system is in a stationary state given by [32] max where we use c = 0.001 and = 1.5 x, and n and m are the two last iteration steps where reinitialization of φ α occurred.
This end state is treated as the quasi-static equilibrium of the system. By incremental changes in boundary conditions and/or constraints the system approximates the quasi-static evolution of a multiphase system.

Two-phase system as a special case
A two-phase system requires only a single level set function φ that evolves according to Eq. (2). In this case we can rescale the velocities from Eqs. (9) and (10) by means of interfacial tension to obtain [31] x S(ψ )|∇ψ| cos β, where β = 180 • − θ . In Eq. (21), C is a prescribed value for capillary pressure divided by interfacial tension, representing interface curvature by Eq. (4). With interfacial tension eliminated from the evolution equation, capillary equilibrium in the pore space occurs when the computed curvature field κ has become equal to the prescribed value C around the interfaces (i.e., C = κ). Note that in the evolution equations of the MLS approach we have chosen to use phase pressures rather than interface curvatures C because we cannot eliminate interfacial tensions from simulations of systems containing more than two fluid phases. Compared to the MLS approach, applying LVC algorithm in a single level set framework requires only a few changes related to the sign of φ. Assuming φ < 0 in phase α and φ > 0 in phase β, we now update interface curvature C a of isolated domains as follows: where it is crucial that s = 1 for a ∈ I α and s = −1 for a ∈ I β . Domain volume and area calculations follow Eqs. (17) and (18), except for a ∈ I β where φ > 0 so that the volume calculation is H(ψ)H(φ)dV .

Parallel aspects of the LVC algorithm
We have implemented the level set methods described in the previous sections within the Structured Adaptive Mesh Refinement Application Infrastructure (SAMRAI) software system [40]. SAMRAI is an object-oriented MPI and C++-based programming framework that provides an extensible and scalable infrastructure for parallel adaptive mesh refinement (AMR) applications [41][42][43].
The implementation strategy for structured AMR employed in SAMRAI is the so-called patch-based approach, as opposed to the quad/octree approach which is the other typical alternative. In SAMRAI a given processor can be associated with one or several patches. We have implemented AMR in our previous level set methods for two-and three-phase capillarycontrolled displacements in porous media as demonstrated in our recent papers [44,32].
While the patch based parallel approach in SAMRAI can be easily exploited for the level set methods in our previous works, it is more challenging when including the above LVC algorithm. The purpose of this section is thus to explain the implementation of certain important aspects of the parallel LVC algorithm. In the present paper we will only consider uniform grids. The inclusion of AMR introduces additional challenges and is currently work in progress.

Parallel algorithm
We now give the main steps involved in the implementation of our LVC method in the parallel MLS setting. This will closely follow the steps from AlgorithmSeq (Section 3.1), but with the inclusion of additional substeps and comments to explain important extensions needed for parallelism. We will discuss some of these steps in more detail in Section 5.2 as well as in Appendix A. AlgorithmPar 1. Initialize the system: (a) On each patch: Use the grassfire algorithm to identify fluid regions of the conserved phase α with local tags unique to each processor. (b) On each processor: Gather the information in (a) from the corresponding patches and communicate information to the master processor.
(c) On the master processor: Construct unique global tags a ∈ I α (t 0 ) assigned to each isolated fluid region using the algorithm described in Appendix A. Construct the set { a : a ∈ I α (t 0 )} with the help of the information from step (a) above. Communicate relevant information back to all processors and patches. (c) On the master processor: Perform Eq. (13) and communicate relevant information back to all processors and patches. 3. Iterate the multiphase system from t n to t n+1 with Eqs. (3), (9) and (10), using the phase pressures from step 2, and then enforce the projection step [35]. This is done in parallel using the patch-based parallelism in SAMRAI.
4. Repartition the system at time step t n+1 , that is, construct the set Tag regions as boundaries if they belong to the phase connected to inlet or outlet boundaries. In parallel this is done in a similar manner as in step 1 (a)-(c) above.

Extend the partitioning
In parallel this is done in a similar manner as in step

Simple illustration of some aspects of the parallel algorithm
Here, we elaborate on essential parts of the parallel algorithm given above using the following simple two-dimensional example containing two isolated regions (see Fig. 3). Two domains need to be identified in this simple example, and in the top row of the figure we show local labels (tags) as well as the various patches for one, eight and sixteen processors, respectively. Obviously, local and global labels coincide in the case of a single processor. Note that in parallel computing local tags unique to each processor are identified by performing the grassfire algorithm on each patch. After communicating the various results to the master processor global tags have to be associated with each isolated region, since a given region can cover several patches (and processors). This is done by using the algorithm described in Appendix A.
In the middle row of Fig. 3 the regions are denoted by 1 and 2 , with corresponding global labels 1 and 2, respectively. Outside of the two domains in these three subfigures we show the distance map, which is used for computing the extended regions. This is just the negative integer distance computed with respect to the nearest domain (i.e. −d( x, )). In patches without domains this negative integer distance is simply put to the negative sum of the number of grid points in each coordinate direction. The corresponding extended regions are shown in Fig. 3, bottom row. In the parallel computing process we iterate through the patches connected to a processor and perform the "extension algorithm" (described in Section 3) for each patch. However, when this is completed we must yet again iterate through the patches and possibly correct the results from the extension based on patch neighbourhood information. In essence we compare the "negative distance map" of neighbouring patches and give preference to the shortest negative distance when updating the "extension map".
As seen in Section 3, volume and area computations of the different regions are central to the LVC algorithm. We perform such computations on each patch in parallel. Different contributions are then gathered on each processor and communicated to the master processor where they can be put together with the help of the global domain identification tags discussed above.
The volume redistribution algorithm described in the previous section is also at first done on each patch in parallel in order to find the various domain sizes and corresponding intersections ("hits"), both represented as integers, between old and new domains. Communications are again needed to gather the various data structures on the master processor. After a somewhat involved process involving matrix-like structures the new domain volumes can be computed as in the sequential algorithm (see Section 3). Finally, the pressures in the different regions can be computed, and then communicated to the various processors for use in the level set algorithm.

Parallel scalability of the LVC algorithm
The parallel version of the LVC algorithm introduces additional MPI communication in the LS/MLS codes. In this section we thus investigate the parallel performance of the LVC algorithm using the MLS method. We present tests for both weak and strong scalability of the parallel code based on simulations of curvature-driven evolution of isolated domains toward a steady state where the domains attain spherical shapes. In most realistic applications the size of the data structures involved in the MPI communication will be relatively modest, but it will grow with the number of isolated domains that needs to be identified. Since the parallel scalability thus clearly will depend on this number, we show results where the number of isolated domains ranges from 0 to 216, representing typical simulations in realistic pore geometries.
In the scaling tests the distributions of the phases subjected to LVC (that is, the oil and gas phases) constitute fixed numbers of isolated cubic domains, whereas the remaining phase (water) is continuous. We consider two such cases, as shown in Fig. 4(a). The first case contains 64 (= 4 3 ) domains (32 of each phase) and the second case contains 216 (= 6 3 ) domains (108 of each phase). We explore scaling behaviour in each of these cases for the settings without LVC and with LVC applied to one (i.e., oil) and two (i.e., oil and gas) phases. To maintain control of the number of isolated domains present in the system the computational domain consists entirely of pore space. In the weak scaling tests we first simulate a computational domain with 40 3 grid cells using a single processor. We then refine this configuration three times (keeping the number of isolated domains constant) so that the problem size and number of processors used in the simulations both increase with a factor of 8 for each refinement. We thus end up with a problem containing 320 3 grid cells in the case of 512 processors. In the strong scaling tests we perform simulations on the largest problem size (320 3 grid cells) using 128, 256, 512 and 1024 processors. For 512 processors, the simulations are identical to the corresponding simulations from the weak scalability test.
In the weak and strong scaling tests each simulation executes 20 iteration steps with each of the three fluid LS evolution equations, during which three reinitialization solves (each consisting of 10 iterations) occur after the 10th and 20th LS iteration steps, respectively. In addition, at the end of each LS iteration-time step, the simulations also include the projection step [35] and, if applicable, the LVC algorithm. Overall, the three-phase simulations carry out 120 iterations, of which 60 are main LS iterations and 60 are reinitialization iterations. These fractions are representative for a typical usage of our code. All computations were performed on the supercomputer Fram provided by UNINETT Sigma2 -the National Infrastructure for High Performance Computing and Data Storage in Norway.
Results for the weak scaling are presented in Table 1. With a perfect (ideal) parallel speedup the absolute timing for each of these simulations should equal the timing obtained for the corresponding simulation (i.e., with an equal number of conserved phases and isolated domains) on a single processor. In the case without LVC we observe a very good parallel performance. However, in our opinion the results are also satisfactory when we include LVC on one or two phases, even for  The results from the strong scaling tests are presented in Table 2 and Fig. 4(b). We observe that the parallel speedup is very good for the case without LVC. For the cases with LVC the parallel speedup decreases with an increasing number of isolated domains present in the system. This is due to parallel communication introduced by the LVC algorithm. However, the parallel speedup obtained is still absolutely acceptable and definitely demonstrates that our codes run efficiently on modern supercomputers.
We demonstrate the LVC algorithm on different pore geometries using a fluid system of gas (g), oil (o) and water (w) with physically realistic interfacial tensions σ g w = 0.03, σ ow = 0.02 and σ go = 0.011 N/m. The characteristic surface tension γ * is 0.01 N/m, which together with Eq. (6) implies that γ w = 1.95, γ o = 0.05 and γ g = 1.05. We will also ensure that all nonzero γ α are sufficiently larger than x, to obtain correct behaviour in the limit γ α → 0. For a discussion of vanishing viscosity solutions with level set methods we refer to [36].
For a three-phase fluid system in thermodynamic equilibrium with a solid the three contact angles θ αβ cannot be specified independently. Instead, Eq. (5) (or, equivalently Eq. (7)) implies that the interfacial tensions and contact angles satisfy σ g w cos θ g w = σ go cos θ go + σ ow cos θ ow . (23) In this work, we specify θ ow and calculate θ go and θ g w from linear cos θ go (cos θ ow )-and cos θ g w (cos θ ow )-relationships [45] that satisfy Eq. (23). Our model uses a corresponding set of intersection angles β α that are consistent with these relationships [32]:

2-D sinusoidal pore
We perform simulations on an idealized 2-D pore throat geometry to investigate effects of grid resolution on conservation of volume (area in 2-D) and pressure convergence of isolated domains. The characteristic quantities are L * = 1 × 10 −4 m and p * = γ * /L * = 100 N/m 2 . In dimensionless variables, we consider the sinusoidal pore throat We use initial gas saturation S g,1 = 0.05 and saturation step S g = 0.02. The converged fluid configuration from the last saturation step is the initial configuration for the next. Note that the saturation-controlled approach mimics quasi-static invasion at constant low flow rate, while it predicts capillary pressure at each stationary state [32,33]. Pressure-controlled displacement, which is the other alternative, predicts saturations for each stationary state given a series of fixed capillary pressures. Fluids can enter and leave the system at the inlet and outlet boundaries (that is, x = 0 or x = 2). We perform simulations with We evaluate domain pressures, p a , a ∈ I α , by means of their pressure difference to the other two continuous phases (i.e., local capillary pressure), which represents interface curvature when scaled by interfacial tension, see Eq. (4).

Double displacement with domain fragmentation and coalescence
First we investigate LVC of oil only, while water pressure is constant. Hence, we assume sub-resolution wetting films along the pore walls connect the water to inlet/outlet boundaries. For this purpose we use θ ow = 20 • . Calculation of the other contact angles based on Eqs. (7) and (24) yields θ go = 24.2 • and θ g w = 16.1 • . We wish to investigate double displacements where a continuous invading phase (e.g., gas) displaces a disconnected phase cluster (e.g., oil) that in turn displaces a continuous third phase (e.g., water). Pore-scale experiments have shown that this mechanism can mobilize oil in capillary-dominated three-phase flow [4,5,15,16]. In our simulation example, gas displaces an oil drop that coalesces with static oil drops ahead of the frontal oil/water interface during the motion, while oil-drop fragmentation forms trapped oil layers in pore corners behind the trailing gas/oil interface. Fig. 5 shows only small grid resolution effects on the stationary fluid configurations from these gas-invasion simulations, except for oil drop coalescence events that in some cases occur at slightly different gas saturations. As expected, the capillary pressures P aw , a ∈ I o , between static oil drops and surrounding water are constant and display excellent convergence behaviour with decreasing x, see Fig. 6(a). For the oil drop displaced by gas, the local capillary pressures P aw and  P ga , a ∈ I o , vary during displacement, as the frontal oil/water and trailing gas/oil interfaces move through wide and narrow pore sections, see Figs. 6(b, c). We obtain excellent capillary pressure convergence in this case as well, despite coalescence and splitting events occur during the displacement. Further, Fig. 6(c), shows that the capillary pressures P ga , a ∈ I o for the oil layers trapped by gas in the pore corners behind the front stabilize at constant values. The capillary pressure between the continuous gas and water phases P g w obtains local maxima when either the frontal oil/water interface or trailing gas/oil interface invades the narrowest pore sections. Although P g w varies with gas saturation, Fig. 6(d) shows that the effects of x on P g w are small.

Figs. 7(a, b) show that the domain volume relative errors E[V a
], a ∈ I o , decrease consistently with decreasing x. For the static oil drops surrounded by water E[V a ] remains constant during gas invasion, see Fig. 7(a), while for the oil drops contacting gas E[V a ] varies more with gas saturation due to displacement, coalescence and splitting events ( Fig. 7(b)). Fig. 7(c) shows that the relative error for the total disconnected oil volume E[V o ] behaves similarly. Examining the behaviour at three gas saturations shows that E[V o ] decreases quadratically with decreasing x, see Fig. 7(d).
We now investigate the above example with LVC applied to both oil and water in the simulations. The sequence of stable fluid configurations presented in Fig. 8 shows that the oil drop displaced by gas still merges with the static oil drops surrounded by continuous water during its motion. However, this time isolated water domains break up from the continuous water and trap in pore corners surrounded by gas or oil behind the frontal oil/water interface. For the lowest grid resolution ( x = 0.02) two oil drops completely isolate some water from its continuous phase, leading to a different displacement pattern than for x ≤ 0.01. Formation of oil layers between bulk gas and corner water is unlikely to occur in this example because the fluid system represents a non-spreading oil (σ g w − σ go − σ ow < 0) with unfavourable contact angles (θ go > θ ow ) and because the amount of water trapped in pore corners is large. Instead, triple junctions form as the trailing gas/oil interface passes water domains trapped in the corners, see Fig. 8. Fig. 9 investigates accuracy of capillary pressure between isolated oil drops a ∈ I o and continuous gas and water phases as x decreases. The set of oil drops not in contact with gas converge to constant capillary pressure levels as before, see Fig. 9(a). However, Figs. 9(b, c) show that the capillary pressures for oil drops contacting gas vary slightly more with x than in the previous case where LVC was applied to oil only. This is also the case for the capillary pressure between continuous gas and water, P g w , as shown in Fig. 9(d). We anticipate this is an effect of preserving the volumes of coarsely resolved water domains trapped in the pore corners, which in turn will affect the pressure calculations of both the water and oil domains as well as the continuous gas. However, the relative errors for the individual oil-drop volumes, E[V a ], a ∈ I o , and the relative error for the total disconnected oil volume, E[V o ], show the same consistent trends as before, see Fig. 10. Fig. 11 presents capillary pressure and volume preservation properties for the disconnected water. The capillary pressures between continuous gas and water domains a ∈ I w exhibit some variations with x due to grid-resolution effects of the preserved water layers in the pore corners, see Fig. 11(a). Further, the relative error of individual water volumes  E[V a ], a ∈ I w , decreases consistently with x, see Fig. 11(b). However, Fig. 11(c) shows that the relative errors for the total disconnected water volume, E[V w ], are similar for x = 0.02 and 0.01. This is because the configurations of water differ significantly; for x = 0.02 a large water domain spans the pore cross-section during a significant part of the gas invasion, while for x ≤ 0.01 all isolated water domains exist as layers in the pore corners. Thus, for x ≤ 0.01, E[V w ] decreases quadratically with decreasing x, as shown in Fig. 11(d).

Displacement chains with multiple domains
Here we consider gas invasion in the sinusoidal pore with LVC of oil and water under a more weakly water-wet state, represented by the consistent contact angle set θ ow = 60 • , θ go = 21.3 • and θ g w = 47.5 • . A large oil/water contact angle (θ ow > θ go ) favours formation of oil layers between bulk gas and water in pore corners and allows the oil drops to contact pore walls and isolate intermediate water domains. Then, invasion of continuous gas can displace a chain of several isolated water and oil domains that in turn displaces the water phase connected with outlet, while oil and water layers deposit in pore corners behind the trailing gas/oil interface. Such displacement chains can mobilize disconnected oil ganglia over multiple water-alternate-gas invasion cycles, as observed in micromodel pore-scale experiments [9][10][11].
In our example (for x ≤ 0.01), Fig. 12 shows that the two leftmost oil drops contact the pore walls and isolate a water domain from the continuous water phase, resulting in the displacement chain g → o → w → o → w. The three rightmost oil drops are static as they are completely surrounded by continuous water. However, during gas invasion, these oil drops coalesce with the nearest oil drop approaching from the left, while the leftmost oil drop and the water domain split to form oil and water layers in pore corners behind the trailing gas/oil interface. Note that the mobile oil drop surrounded by isolated water to the left and connected water to the right has different interface curvature (that vary during displacement) on the left and right sides due to significantly different pressures in the water domain and the continuous water phase.
Further, Fig. 12 shows that simulations with x ≤ 0.005 capture configurations with oil layers sandwiched between corner water and bulk gas, whereas the coarse-grid simulation ( x = 0.02) exhibits significant deviations to the other results. Fig. 13 explores accuracy of the capillary pressures between disconnected oil drops a ∈ I o and continuous gas and water phases with decreasing x. Fig. 13(a) shows that oil/water capillary pressures for the subset of oil drops surrounded by water converge as x decreases for both the static oil drops where capillary pressures are constant and the mobile oil drops where capillary pressures vary with displacement. Figs. 13(b) and (c) demonstrate convergence of capillary pressures with decreasing x for the oil drops contacting gas. However, a few scattered data points deviate from this behaviour due to coarse grid resolution of oil layers in pore corners. For x = 0.005 oil layers break into two oil drops on the gas/water interfaces as the trailing gas/oil interface moves into narrower pore cross-sections accompanied by higher gas pressure. These oil drops attain a high and negative gas/oil capillary pressure and a corresponding high and positive oil/water capillary pressure. The oil layers repair themselves with normalized domain capillary pressures when the trailing gas/oil interface moves into wider pore regions accompanied by lower gas pressure. For x = 0.0025 oil layers are unaffected by these gas pressure variations that depend on the position of the trailing gas/oil interface. Fig. 13(d) shows that the capillary pressure between continuous gas and water phases, P g w , converges consistently with decreasing x.    Fig. 15 presents the results for the water domains. The capillary pressures between continuous gas and the water domains converge with decreasing x. However, a few small water layers with low grid resolution give a couple of deviating data points. The water domains exhibit a similar volume-conservation behaviour with decreasing x as observed for the oil phase.

2-D synthetic porous medium
Here, we test the LVC algorithm on a 2-D synthetic porous medium with characteristic length L * = 1 × 10 −4 m, consisting of 110 circles (representing solid phase) with dimensionless radii in the range [0.04, 0.07]. The circle centres are placed at uniformly distributed random positions within the computational domain ∈ [0, 1] × [0, 2] without eliminating possible circle overlaps. The space outside the circles represents pore space. First we use the two-phase LS method with LVC of water and simulate oil invasion from y = 0 into a fully water-filled pore space (primary drainage), for both θ ow = 20 • and 60 • . From three different oil saturations established during primary drainage we simulate gas invasions at constant P ow , also from y = 0, using MLS method with LVC of both oil and water. For primary drainage, we simulate saturation-controlled oil invasion with S o = 0.02 and update oil/water interface curvature continuously according to given a set of target oil volumes V target o,k that correspond to target oil saturations S o,k . The following gas invasions occur in a similar manner, using S g = 0.01 and Eq. (27). For the LS functions, we use reflective boundary conditions except at the inlet ( y = 0) where we add a pore-space layer with thickness 2 x occupied by invading phase and use linear extrapolation.
For the phases subjected to LVC, we assign the continuous phase pressure to isolated domains connected to boundaries y = 0 (inlet) and y = 2 (outlet), while we preserve domains contacting the other boundaries (x = 0 and x = 2). The LS functions are reinitialized, according to Eq. (19), for each 10th iteration steps.   Fig. 12 (d). The left columns of Figs. 16 and 17 show the initial oil and water configurations prior to gas invasion for θ ow = 20 • and 60 • , respectively. As we expect, the amount of water left behind the front during primary drainage is higher for θ ow = 20 • than for θ ow = 60 • although both simulations trap a large water domain in the lower left part of the pore geometry.
These configurations also show that the oil phase has a different interface curvature to the continuous water phase than to the various disconnected water regions. Fig. 18 shows the distribution of domain capillary pressures compared with the continuous-phase capillary pressure curve (in physical units). The capillary pressures for isolated water domains are significantly lower for θ ow = 60 • than for θ ow = 20 • . Several water domains stabilize at negative capillary pressures due to the higher contact angle.
Gas invasion from the two lowest oil saturations occurs by double displacement of a large disconnected oil cluster, while smaller oil domains and layers get trapped between disconnected water and continuous gas, see Figs. 16(a, b) and 17(a, b). Branching of the gas phase triggers another double displacement of oil toward outlet boundary at the upper right part of the pore geometry. The number of trapped oil layers left behind the invading gas is higher for θ ow = 60 • than for θ ow = 20 • because oil invades narrower pore constrictions during primary drainage under weak water-wet states. Gas invasion at high initial oil saturation occurs by direct gas/oil displacement, while fragmentation of continuous oil traps oil domains behind the front, see Figs. 16(c) and 17(c). The extent of displacement chains with multiple domains is limited in these simulations. The most conspicuous example is displacement g → o → w → o, as highlighted by an arrow in the fourth image of Fig. 17(c). This is consistent with micromodel studies where such multiple displacements occur more frequently for subsequent water-and gas-invasion cycles [9][10][11].
The three-phase configurations in Figs. 16  has high gas/water contact angle (θ g w = 47.5 • ), more water domains have lower (and even negative) capillary pressure than the case with θ ow = 20 • where θ g w = 16.1 • . Figs. 19(c) and (d) present gas/oil capillary pressure with respect to continuous oil and disconnected oil domains for the gas invasion simulations from endpoint oil saturations 0.63 (θ ow = 20 • ) and 0.72 (θ ow = 60 • ) after primary drainage. Note that several domain capillary pressures are negative despite gas/oil contact angles are small (θ go = 21.3 • and 24.2 • ). This is possible when gas/oil/water triple lines with water cusps, rather than gas/oil/solid contact lines, form on the pore walls. For such instances, configurations of isolated oil domains resemble oil drops (confined by pore geometry) resting on gas/water interfaces with negative gas/oil capillary pressures and positive oil/water capillary pressures. The ranges of gas/oil capillary pressures for oil domains are similar for both θ ow = 20 • and θ ow = 60 • because the gas/oil contact angles in the two cases are about the same.

3-D porous medium
We use a pore structure extracted from a segmented microtomography data set of Castlegate sandstone [46] to explore effects of LVC in simulations of capillary pressure curves on realistic 3-D porous media. The subset contains 128 3 voxels, where physical voxel length is 5.6 μm. Thus, L * = 7.168 × 10 −4 m. The simulations use contact angles θ ow = 40 • , θ go = 23.1 • , and θ g w = 32 • . We use reflective boundary conditions for the LS functions at all faces of the computational domain, except at inlet (z = 0), where we add a pore-space layer (thickness 2 x) occupied by invading phase and use linear extrapolation. For phases with LVC we assign the continuous-phase pressure to domains connected with top or bottom boundaries of the sample, while we preserve domains contacting side walls. We simulate pressure-controlled displacement  First we investigate effects of LVC applied to the wetting phase (water) during primary drainage, using the two-phase LS method (Section 4). Fig. 20(a) shows that LVC has little effect on capillary pressure in the high water saturation regime, but substantial differences occur for S w < 0.5. This leads to significantly different initial water saturations S wi (that is, endpoint saturations) in the two simulations. We obtain S wi = 0.23 with LVC in contrast to S wi = 0.02 without LVC. Figs. 20(b) and (c) show that the difference occurs because large water domains form midway in the z-direction in the case with LVC. Previously, we performed LS simulations without LVC on the same sample to investigate the impact of using one level with adaptive mesh refinement (AMR) around fluid interfaces [44]. AMR raised the P ow -level in the low S w -range, but S wi remained about the same. Simulations with θ ow = 0 • , 20 • , and 40 • , gave S wi = 0.05, 0.04, and 0.02, respectively [44].
In comparison, core-scale measurements of S wi on homogeneous sandstone typically vary between 0.02 and 0.13 [47,48].
Specifically, Shikhov et al. [48] measured S wi = 0.11 on Castlegate sandstone. Apart from the difference in rock volumes, capillary pressure curves from pore-scale simulations and core-scale measurements typically deviate at low wetting-phase saturations because of a lacking representation of microporosity on the rock   image and because simulation models do not capture thin sub-grid-scale wetting-phase films along the pore walls [49][50][51][52][53][54]. In our simulations, the case without LVC assumes all water is continuous through thin films, whereas the case with LVC traps all water domains that appear isolated. We anticipate AMR simulations combined with higher image resolution can allow for drainage to a lower S wi with LVC because the existence of more water filaments or layers in pore corners connects water domains and makes water more continuous. However, incorporating AMR within the current LVC algorithm requires further development. Another aspect is that calculations of S wi might depend on whether phase domains contacting side walls shall preserve their volumes or not. We expect that these LVC boundary conditions become less sensitive to S wi with larger computational domain due to a higher volume/area ratio.
Finally, we compare simulations with LVC applied to oil against simulations with LVC applied to both oil and water for gas invasions from S w = 0.67 after primary drainage. The simulations end when the oil saturation no longer decreases. Fig. 21 displays minor differences between the saturation paths, whereas the capillary pressure between continuous gas and water is consistently higher for the case with LVC applied to both oil and water. The case with LVC of oil reaches a lower oil saturation, suggesting that volume conservation of water makes it harder to displace disconnected oil drops. Fig. 22 shows three-phase fluid configurations during the gas invasions. In both simulations, gas displaces continuous oil directly and disconnected oil by double displacement. These mechanisms also lead to oil fragmentation with trapping of smaller oil drops in narrower pore spaces.

Summary and conclusions
This work presents a multiphase level set method with LVC for capillary-controlled pore-scale displacements in porous media. The method conserves volumes of disconnected phase ganglia and computes capillary pressures between ganglia and surrounding phases. A conservative volume redistribution technique handles ganglia splitting and coalescence. The algorithm is extended to handle domains that spanned multiple patches for simulation using parallel codes.
Although the method allows for an arbitrary number of phases, here we perform simulations of two-and three-phase displacements with LVC of one and two phases. Specifically, we provide simulation examples of important three-phase displacement mechanisms in porous media, like double displacements and displacement chains of multiple domains with splitting and coalescence.
Simulations on idealized 2-D pore geometries with increasing grid resolution demonstrate convergence of ganglia capillary pressures as well as consistently decreasing relative error of individual ganglia volumes. Relative error for the global disconnected phase volume decreases quadratically with decreasing grid-cell length. Simulations of gas invasion after primary drainage in 2-D water-wet porous media show that double displacement with oil fragmentation is a dominant mechanism. Disconnected oil and water ganglia develop large capillary pressure ranges that encase the simulated capillary pressure between corresponding continuous phases. Static oil drops surrounded by continuous water have constant capillary pressures, while the capillary pressures of oil drops contacting gas and water fluctuate during double displacement through wide and narrow pore regions. Finally, we demonstrate that the main effects of using LVC in level set simulations of primary drainage and gas invasion on 3-D porous rock is higher capillary pressure level for the continuous phases and higher endpoint saturations.
Future work will use the developed method to investigate residual phase trapping as well as pressure states and topology of isolated oil and gas ganglia during gas and water invasion cycles for recovery and storage processes in porous rock. We also anticipate that using the local conservative level set approach for interface tracking as an alternative to standard level set methods has potential to improve mass conservation significantly in multiphase flow simulations based on the Navier-Stokes equations. In the LVC algorithm we compute global labels to identify the various domains, see algorithm in Section 5. If we were to perform the grassfire algorithm in one go it would demand a lot of communication between processors, and this could be a potential bottleneck for the computational efficiency. Instead we perform the grassfire algorithm on each individual patch connected to a given processor. We make sure that each processor will a unique labelling of its distinct regions. We call these regions for subdomains as it is possible that they can be part of domains that span multiple processors.

CRediT authorship contribution statement
The purpose of the algorithm we describe here is to identify which subdomains that are part of the same domain. To do this we need unique tags for the subdomains, and that each processor makes a "list of neighbours" that contains which subdomains are in contact through the boundaries between the processor and its neighbours. Each processor sends its "list of neighbours" to a master processor where the algorithm for identifying connected subdomains is preformed.
We use the ordered pair, a = (a r , a l ), of the processor rank, a r and the local label, a l , as a subdomain's unique tag. If subdomains with tags a and b are connected then we will add the pair [a, b] to a "list of neighbours". As an example, the "list of neighbours" for processor with rank 3 in where pair of connected subdomains are in square parenthesis. In the algorithm we will also need an ordering of the tags, which we define so that a < b if a r < b r , or if a r = b r and a l < b l further, a = b if a r = b r and a l = b l .
Since all subdomain tags are represented as the 1st element in the pairs in L, we can use it to define a function, f L , so given that (a, b) ∈ L then f L (a) = b, e.g. f L ((2, 1)) = (2, 1) for the initial L in the example above. Note that f L is redefined each time L is changed. Our goal is now to identify the set of all subdomains that either directly or indirectly share a boundary. To do this we will define a function g L so that, when all "list of neighbours" are added to L, all subdomains part of the same domain have the same value of g L or more precisely g −1 L (g L (a)), where g −1 L is the inverse of g L , will be the set of all the tags of subdomains that is part of the same domain as a. To accomplish this we use f L to define a recurrence relation a (n+1) = f L (a (n) ), where a = a (0) , and we will update L so that the sequence f L generates is decreasing, that is a (n+1) ≤ a (n) . We remark that n ≥ 0 and limited by the total number of subdomains. For instance, if the sequence only contains one element we have simply identified an isolated subdomain. We note that this is fulfilled by the initial L, since a = f L (a) for all a ∈ L. The sequence a (n) is bounded and will attend its minimum value which is a fixed point, a * , of f L . We define the function g L so that g L (a) = a * .
For each entry, [a 1 , a 2 ], in the "lists of neighbours", sent to the master processor, we have the values a * 1 = g L (a 1 ) and a * 2 = g L (a 2 ). If a * 1 < a * 2 , we change the tuple (a * 2 , a * 2 ) in L to (a * 2 , a * 1 ), so that f L (a * 2 ) = a * 1 , if a * 1 > a * 2 , we change L so that f L (a * 1 ) = a * 2 , or if they are equal L is unchanged. We note that this redefinition of L is in accordance with the assumption of a decreasing sequence generated by the recurrence relation. When this is done for all "lists of neighbours", the function g L partition the subdomains into domains, where subdomains with the same value of g L is part of the same domain. We can justify this by noting that when a pair of subdomains from a "list of neighbours" with tags a 1 and a 2 are added to L, L is changed so that g L (a 1 ) = g L (a 2 ) = a * , and we have the sequences, a 1 > a (1) 1 > . . . , a (n 1 ) 1 = a * and a 2 > a (1) 2 > . . . , a (n 2 ) 2 = a * . When a new pair of subdomains, from a "list of neighbours", is added, we change L to L . If the pair contains tags equal to elements in any of the sequences a (n) 1 and a (n) 2 defined above, we know that it is only the value of f L (a * ) that can change, per definition, so we still have that g L (a 1 ) = f L (g L (a 1 )) = f L (g L (a 2 )) = g L (a 2 ), and by iterating this argument we know that they will have the same value when all pair of tags in all "lists of neighbours" are added. So as long as two subdomains are related through a sequence of subdomains they will have the same value of g L and thus be part of the same domain.
We have not exhausted all possible pairs in "the lists of neighbours", but we note that the domain structure is recognized. For the first four processors we have identified a domain that spans three of processors.