A Monte Carlo Algorithm for Immiscible Two-Phase Flow in Porous Media

We present a Markov Chain Monte Carlo algorithm based on the Metropolis algorithm for simulation of the flow of two immiscible fluids in a porous medium under macroscopic steady-state conditions using a dynamical pore network model that tracks the motion of the fluid interfaces. The Monte Carlo algorithm is based on the configuration probability, where a configuration is defined by the positions of all fluid interfaces. We show that the configuration probability is proportional to the inverse of the flow rate. Using a two-dimensional network, advancing the interfaces using time integration the computational time scales as the linear system size to the fourth power, whereas the Monte Carlo computational time scales as the linear size to the second power. We discuss the strengths and the weaknesses of the algorithm.


Introduction
The characterization of porous media at the pore level is undergoing a revolution [1]. Through the use of new scanning techniques, we are capable of reconstructing the pore space completely, including the tracking of motion of immiscible fluids. A gap is now appearing between the geometrical characterization of porous media and our ability to predict their flow properties based on this knowledge.
The pore scale may be of the order of microns whereas the largest scalese.g. the reservoir scale -may be measured in kilometers. Hence, there are some eight orders of magnitude between the smallest and the largest scales. At some intermediate scale, that of the representative elementary volume (REV), the porous medium may be regarded as a continuum and the equations governing the flow properties are differential equations. The crucial problem is to construct these effective differential equations from the physics at the pore scale. This is the upscaling problem. A possible path towards this goal is to use brute computational power to link the pore scale physics to pore networks large enough so that a continuum description makes sense. Alas, this is still beyond what can be done numerically. However, computational hardware and algorithms are steadily being improved and we are moving towards this goal.
It is the aim of this paper to introduce a new algorithm that improves significantly on the efficiency of network models [2]. These are models that are based on the skeletonization of the spaces in such a way that a network of links and nodes emerge. Each link and node are associated with parameters that reflect the geometry of the pore space they represent. The fluids are then advanced by time stepping some simplified version of equations of motion of fluid. The bottle neck in this approach is the necessity to solve the Kirchhoff equations to determine the pressure field whose gradients drive the fluids in competitions with the capillary forces.
A different and at present popular computational approach, among several, is the lattice Boltzmann method [3,4]. This method, based on simultaneously solving the Boltzmann equations for different species of lattice gases, is very efficient compared to the network approach necessitating solving the Kirchhoff equations. However, the drawback of the lattice Boltzmann approach is that one needs to resolve the pore space. Hence, one needs to use a grid with a finer mask than the network used in the network approach. This makes the lattice Boltzmann approach very efficient at the scale where the actual shape of the pores matter, but not at the larger scale where the large scale topology of the pore network is more important. Further methods which resolve the flow at the pore level are e.g. smoothed particle hydrodynamics [5,6,7] and density functional hydrodynamics [8]. When network models are so heavy numerically that the networks that can be studied are not much larger than those studied with the pore scale methods, the latter win as they can give a more detailed description of the flow. However, if the computational limitations inherent to network models could be overcome, they would form an important tool in resolving the scale-up problem: at small scale network models would be calibrated against the methods that are capable of resolving the flow at the pore level. On large scales, their results may be extrapolated to scales large enough for homogenization, i.e., replacing the original pore network by a continuum.
As pointed out above, the bottleneck in the network models is the necessity to determine the pressure field at each time step. When the time steps are determined by the motion of the fluid interfaces, these will be small as they typically are set by the time lapse before the next interface reaches a node in the network. Time stepping allows detailed questions concerning how flow patterns develop in time to be answered. That is, the time stepping provides a detailed sequence of configurations where each member of the sequence is the child of the one before and the parent of the one after. If the quantities that are calculated are averages over configurations, time stepping will provide too much information; for averages the order in which the configurations occur is of no consequence. If the order in which the fluid configurations occur is scrambled, the averages remain unchanged. This is where the Monte Carlo method enters. It provides a way to produce configurations that will result in the same averages as those obtained through time stepping. The order in which the configurations occur will be different from those obtained by time stepping. The time stepping procedure necessitates that there are tiny differences between each configuration in the sequence, since the time steps have to be small. This limitation is overcome in the Monte Carlo method which we will describe here. This makes the Monte Carlo method much more efficient than time stepping as we will see.
In Section 2 we describe the network model we use to compare the Monte Carlo method with time stepping, see Aker et al. and Knudsen et al. [9,10]. In the next Section 3, we start by explaining the statistical mechanics approach to immiscible two-phase flow in porous media that lies behind the Monte Carlo algorithm we propose [11,12]. In particular, we derive the configuration probability -the probability that a given distribution of fluid interfaces in the model will appear. This is also known as the ensemble distribution in the statistical physics community. Based on this knowledge, we then go on to describe the Monte Carlo algorithm itself. This section is followed by Section 4 where we compare the Monte Carlo method with time stepping using the same network model described in Section 2. We then go on to compare the efficiency in terms of computational cost of the two methods. We end this section by discussing the limitations of the Monte Carlo algorithm as it now stands and point towards how these may be overcome. We end by Section 5 where we summarize the work and draw our conclusions.

Network Model
In order to have a concrete system to work with, we describe here the details of the network model we use. The model is essentially the one first developed in references [9,10]. For simplicity we do not consider a reconstructed pore network based on a real porous medium [13,14]. Rather, we simply use a two-dimensional square network, with disorder in the pore radii, oriented at 45 • with respect to the average flow direction as shown in Figure 1. As described in [10], we use bi-periodic boundary conditions. Hence, the network takes a form of the surface of a torus. In this way, the two-phase flow enters a steady state after an initial transient period. This steady state does not mean that the fluid interfaces are static. Rather, we use capillary numbers high enough so that fluid clusters incessantly form and break up. By steady state we mean that the macroscopic averages -averages over the entire network -are well defined and do not drift.
The network contains L × L links. All links have equal length l, but their radii have been drawn from a uniform distribution of random numbers in the interval [0.1l, 0.4l]. We set l = 1mm. We neglect gravitational effects.
Fluid flow through each link in the network is modeled using the Washburn equation [15], see Figure 2. There is a volume flow q passing through it driven by the two pressures p 1 and p 2 . Each fluid interface contributes a capillary pressure p c (x) where x ∈ [0, l] is the position of the interface. The capillary pressure is given by the Young-Laplace equation where γ is the surface tension, θ the contact angle between the interface and the pore wall. We set γ cos θ = 30 dyn/cm. r 0 is the average link radius. We assume that the link has a shape so that p c attains the given x dependence. It has been chosen so that p c (0) = p c (l) = 0 and max x |p c (x)| = |p c (l/2)|. The Washburn equation then becomes where µ av = s nw µ nw + s w µ w is the viscosity. s nw = l nw /l and s w = l w /l are the fractions of the link length that cover the non-wetting and wetting fluids respectively so that s nw + s w = 1. We set µ nw = µ w = 1 poise. We define the capillary number Ca as where · · · is an average over all links. A pressure difference ∆P is applied across the network. This is done in spite of the network being periodic in the direction of the pressure difference, see Knudsen et al. [10]. By demanding balance of flow at each node using the Washburn equation (2), we determine the pressures (p i ) at the nodes. This is done by solving the q l x 0 p p 1 2 p(x) c corresponding matrix inversion problem by using the conjugate gradient algorithm [16].
When the pressures at nodes are known, the flow q ij -here between neighboring nodes i and j connected by a link -is calculated using equation (2). Knowing the velocity of the interfaces in each link, we then determine the time step such that any meniscus can move a maximum distance, say, one-tenth of the length of corresponding link in that time. All the interfaces are then moved accordingly and the pressure at the nodes are determined again by conjugate gradient algorithm. This is equivalent to event-driven molecular dynamics. When an interface reaches the node, the interface will spread into the links that are connected to the node and which have fluid entering them from the node. The rules for how this is done are described in detail in Knudsen et al. [10].

Metropolis Monte Carlo
We first describe the theory that lies behind the Monte Carlo algorithm that we present. We need to introduce the concepts of configuration, and configuration probability, also known as the ensemble distribution in the statistical mechanics community. We then go on to derive the configurational probability. Armed with this, we construct the Monte Carlo algorithm [17] after having presented a short review of the Metropolis version of Monte Carlo [17,18].

Statistical Mechanics of Immiscible Two-Phase Flow
Sinha et al. [19] studied the motion of bubbles in a single capillary tube with varying radius. Suppose that the capillary tube has a length L and a radius that varies as r = r 0 /[1−a cos 2πx/l] where l ≪ L, a is an amplitude and r 0 the average radius of the tube. Suppose furthermore that the tube is filled with wetting fluid except for a bubble of length ∆x b and a center position x b . By using equation (1), one derives the net capillary force from the two interfaces that limit the bubble as, where σ = 4aγ cos θ sin(π∆x b /l)/r 0 . By combining this equation with the Washburn equation (2), one findṡ where πr 2 0ẋb = q, and ∆p = (L/l)∆P where ∆P is the pressure difference across the tube.
Suppose there is a quantity f = f (x b ) that depends on the position of the bubble in the capillary tube. For example, f might be the flow q. Let us now assume that ∆P does not vary in time. The time average of f is then where x b (t) is the time integration of the Washburn equation (2) and the time period We note, and this is the crucial observation, that we may change integration variable from time t to bubble position where is the configuration probability. That is, the configuration of the tube is given by the position of x b of the bubble. Equation (8) gives the probability density to find the bubble at position x b -and hence in that configuration. The Washburn equation (5) gives the motion of the bubble that is used in equations (7) and (8). The Washburn equation assumes that we control the pressure drop ∆P . If we on the other hand control the flow q, the equation of motion becomesẋ The time period now becomes and hence the configurational probability is which states that all positions of the bubble is equally probable.
To ramp up the complexity of the problem, we assume that there are N bubbles in the one-dimensional tube. The centers of mass of bubble number j ∈ [1, N ] is x j and it has a width of ∆x j . Since the system is one dimensional, all bubbles move with the same speedẋ j =ẋ 1 . The Washburn equation is theṅ where δx j = x j − x 1 and Solving the equations of motion (12) gives x j = x j (t). We may invert , analogous to the one introduced in equation (7). Its time average is where where q = πr 2 0ẋ1 . This is precisely the same expression as in (8). We now turn to complex network topologies. For concreteness, we may imagine a two-dimensional square network. However, the arguments presented in the following are general. A configuration is given by the position of all interfaces. Let us Hence, x i contains information both on which link the interface sits in and where it sits in the link. A flow Q passes through the network. The flow equations for the network consist of a Washburn constitutive equation for each link combined with the Kirchhoff equations distributing the flow between the links. The motion of the interfaces are highly non-linear, but of the formẋ i = g i (x). Solving these equations gives x j = x j (t).
Again we consider a function f = f (x) of the position of the interfaces. Its time average is Here we have inverted x i = x i (t) so that we have t = t(x i ) and then substituted The configurational probability is defined as before, Let us now choose x i = x 1 to be an interface moving in a link that carries all the flow in the network. Such a link is a capillary tube connected in series with the rest of the network. In this case we haveẋ 1 = Q/πr 2 0 , where Q is the total flow. Hence, we have We have in the discussion so far compared the time evolution of a given sample defined by an initial configuration of interfaces. We now imagine an ensemble of initial configurations of interfaces. Each sample evolves in time and there will be a configurational probability (18) for each. This will have the same value for each configuration x that corresponds to the same flow Q. Hence, we have the configurational probability This equation is the major theoretical result of this paper: all configurations corresponding to the same Q are equally probable. Intuitively, equation (19) makes sense: The slower the flow, proportionally the more the system stays in -or close to -a given configuration [20].
Is the system ergodic? Equations (7), (14) and (16) answer this question positively. Time averages give, by construction, the same results as configurational averages.

Implementation of the Metropolis Algorithm
In order to present the details of the Metropolis Monte Carlo algorithm that we propose, we first review the general formulation of the Metropolis algorithm [21,22].

General Considerations
We have a set of configurations characterized by the variable x, the positions of the interfaces. We now wish to construct a biased random walk through these configurations so that the number of times each configuration is visited -i.e., the random walk comes within dx of the configuration -is proportional to Π(x). proportional to the probability for that configuration. The Metropolis algorithm accomplishes this goal. In order to do so, a transitional probability density from state x to state x ′ is constructed as where π(x, x ′ ) is the probability density to pick trial configuration x ′ given that the system is in configuration x. It is crucial that π(x ′ , x) is symmetric, Equations (20) and (21) ensure detailed balance, Detailed balance guarantees that the biased random walk visits the configurations x with a frequency proportional to Π(x). The generated configurations follow the ensemble distribution. When we combine equations (19) and (20), we have

The Implementation
The Metropolis Monte Carlo algorithm based on equation (23) consists of two crucial steps. The first step consists in generating a trial configuration and the second step consists in deciding whether to keep the old configuration or replacing it with the trial configuration. The first step, generating the trial configuration, is governed by the trial configuration probability π(x, x ′ ) which must obey the symmetry (21). That is, if the system is in configuration x, the probability to pick a trial configuration x ′ must be equal to the probability to pick as trial configuration x if the system is in configuration x ′ .
Suppose the system is in configuration x. One needs to define a neighborhood of configurations among which the trial configuration is chosen. If the neighborhood is too restricted, the Monte Carlo random walk will take steps that are too small and hence would be inefficient. If, on the other hand, the neighborhood is too large, the random walk ends up doing huge steps that will miss the details.
We propose generating the trial configurations as follows. Our system is shown in figure 3 and consists of L × L links as described in Section 2. There is a flow q ij through link ij connecting the neighboring nodes i and j. There is a total flow rate Q in the network given by and a corresponding pressure drop ∆P . We choose a randomly positioned sub network as shown in figure 3. The network consists of Λ×Λ links. We "lift" the sub network out of the complete network and fold it into a torus, i.e, implementing bi-periodic boundary conditions. The configurations of fluid interfaces in the sub network remains unchanged at this point.
We calculate the flow rate in the sub network By solving the Kirchhoff equations on the sub network, we time step the configuration forwards in time while keeping the flow rate Θ constant. We end the time integration when 4 -arbitrarily chosen -sub network pore volumes have passed through it. The bi-periodic boundaries of the sub network is then opened up and the sub network with the new configuration of fluid interfaces is placed back into the full network. This is then the trial configuration x ′ .
Part of the probabilistic choice of the trial configuration that defines π(x, x ′ ) rests on the choice of the sub network: its position is picked at random. Hence, if the system is in state x or in trial state x ′ , the probability to pick a particular sub network is the same. This makes this part of the choice of trial configuration symmetric. When the sub network is time stepped for 4 sub system pore volumes, this is done at constant flow rate Θ. Hence, all sub network configurations are equally probable, see equation (19). Hence, also this part of the choice of trial configuration is symmetric. The full probability π(x, x ′ ) is the probability of picking a given sub network times the probability that a given configuration will occur. Combining the two leads to the necessary symmetry (21). We point out here that whereas the configurational probability Π(x) in (19) is valid for all configurations, through the way we generate our samples, we are restricting ourselves to physically realistic samples in that they are generated through time stepping parts of the system. We cannot at this stage prove that this does not bias our sampling.
Once the trial configuration x has been generated, it is necessary to calculate the total flow rate Q = Q(x ′ ) in the network. We then decide to accept the trial configuration x ′ by using (23). This defines a Monte Carlo update.
We repeat this procedure until each link in the network has been part of at least one sub network. This defines a Monte Carlo sweep.

Results
We now present numerical results of the Monte Carlo simulation considering the model described in Section 2 and we will compare them with the results by time stepping simulations. Simulations are performed for two different ensembles, one is when the total flow rate Q is kept constant (CQ ensemble) and the other when the total pressure drop ∆P is kept constant (CP ensemble).
In other words, all trial configurations are accepted.
In figure 4 we plot F nw -the non-wetting fractional flow -as a function of S nw -the non-wetting saturation -where the circles and the squares denote the results from Monte Carlo and time stepping, respectively. The plots, as expected, show an S-shape. This is because the two immiscible fluids do not flow equally, and the one with higher saturation dominates. Hence, the curve does not follow the diagonal dashed line, which corresponds to F nw = S nw , shown in the figure. Rather, F nw is less than S nw for low values of S nw and higher than S nw for higher value of S nw . It therefore crosses the F nw = S nw line at some point, which is not at S nw = 0.5. This is due to the asymmetry between the two fluids, as one is more wetting than the other with respect to the pore walls. This behaviour is more prominent for the lower value of Ca, as capillary forces play a more dominant role. numbers as a function of S nw are shown in figure 5. Similar to the fractional flow plots, we see that the results are same for Monte Carlo and time stepping for a wide range of S nw . We only see differences at high values of S nw . ∆P increases with S nw , reaching a maximum at some intermediate saturation and then decreases again. When S nw increases from zero, more and more interfaces appear in the system causing an increase in capillary barriers associated with interfaces. As the total flow rate Q is constant, a higher pressure is needed to overcome the capillary barriers. The decrease of ∆P after the maximum is due to the decrease of the number of interfaces blocking the fluids.

Constant ∆P ensemble
We now turn to the constant pressure ensemble. Here we keep ∆P constant throughout the calculations. In this case, the Metropolis Monte Carlo algorithm, equation (23), becomes Results for the simulations with constant ∆P are shown in figures 6 and 7. Simulations are performed for two different values of ∆P , 15kPa and 6.5kPa. The steady-state values of F nw show similar variation with S nw as in the constant Q ensemble and we see good agreement between the results for Monte Carlo and time stepping for a wide range of S nw . Here Q varies with the saturation and the corresponding capillary numbers are plotted in figure 7 for Monte Carlo and time stepping. As discussed before, the number of interfaces first increase with the increase in saturation from zero, reaches a maximum value, and then decreases again as S nw approaches 1. The pressure is constant here, so the total flow rate decreases with increasing capillary barriers at the interfaces and correspondingly   Ca varies as in figure 7. Here again, good match between the results Monte Carlo and time stepping can be observed. We show in Table 1 the percentage of rejections for the data shown in Figure  7. The number of rejections is in all cases quite small. This can be understood as follows. Set Q(x, ∆P ) = Q and Q(x ′ , ∆P ) = Q + δ where δ may be positive or negative. Hence, the probability to accept the new configuration is where we have assumed δ ≪ Q. With a small δ the probability to reject the trial configuration is small. This is reflected in Table 1.

Computational Cost
Here we present a detailed comparative analysis of the computational cost of the two algorithms. We do this by measuring the computational time (T MC for the Monte Carlo method and T TS for the time stepping method respectively) for different system sizes L. We use the conjugate gradient method to solve the Kirchhoff equations. This is an iterative solver. When the network contains L × L links (L 2 /2 nodes), each iteration demands L 2 /2 operations. The number of iterations necessary to solve the equations exactly scales as L 2 , making the total cost scale as L β , where β = 2 + 2 = 4. However, in practice, the number of iterations necessary to reach the solution of the Kirchhoff equations to within machine precision is much lower than that needed for the theoretically exact solution. As we shall see, β is much smaller than four.
The number of time steps needed to push one pore volume through the network is n t . We expect it to depend on L as n t = aL τ , where a is a prefactor essentially measuring the number of time steps on the average it takes for an interface to cross a link. In our calculations, this is of the order of 10. Intuitively, this number should be proportional to the width of the network, L, making τ = 1. In practice, as we shall see, it is slightly larger.
For each time step, the conjugate gradient demands t cg = bL β operations where b is another prefactor. The total computational time (T TS ) per pore volume is then where α TS = τ + β. Based on the theoretical considerations above, setting β = 4 and τ = 1, we have T TS ∼ L 5 . The actual computational time measured using Here the time stepping procedure is run for 100 injected pore volumes and in the Monte Carlo method, we do 25 sweeps. Each update is based on running the sub system for 4 injected sub network pore volumes. In this way, when Λ = L = 20, the timing of the two methods are equal. We use the CQ ensemble with Ca = 0.1 and S nw = 0.4. Six different system sizes, L = 20, 40, 60, 80 and 120 are considered. From the slopes, the exponents α for time stepping and Monte Carlo are found. For Monte Carlo, we find α MC = 1.98 ± 0.01, which is close to the theoretically expected value T MC ∼ L 2 (see text). However, for time stepping, we find α TS = 3.99 ± 0.03 which is much smaller than theoretical expectation -T TS ∼ L 5 . In the inset, we plot the average time, t cg , taken by the conjugate gradient solver to solve one entire pressure field. We find t cg ∼ L 2.88±0.02 . The number of time steps per pore volume, n t , scales as n t ∼ L 1.11±0.03 . Combining these two results, we find that the computational time for the time stepping procedure to scales as T TS ∼ L 3.99 .
the clock() function in C is plotted in figure 8 for Ca = 0.1 and S nw = 0.4. We find that T TS scales with L with an exponent α TS = 3.99 ± 0.03 which is much smaller than 5. Measuring n t and t cg independently gives τ = 1.11 ± 0.03 and β = 2.88 ± 0.02, see the insert in figure 8. For the Monte Carlo algorithm, each sweep ideally contains (L/Λ) 2 individual Monte Carlo updates. Each Monte Carlo update consists of time stepping a sub lattice of size Λ×Λ. Hence, the cost of a Monte Carlo update is abΛ α TS when using equation (29). However, each time stepping of a sub lattice is followed by solving the Kirchhoff equations for the entire lattice in order to determine Q for the trial configuration. The cost of this operation is bL β . The time per Monte Carlo sweep is then where α TS = 3.99 and β = 2.88. The factor "4" signifies that we time step the sub lattice for four pore volumes. By setting a ≈ 10 and Λ = 20, the first term will dominate compared to the second term on the right hand side of this equation if 4aΛ α TS ≈ 6.4 × 10 6 > L 2.88 or L > 230 where the second term, which scales as L 4.88 , starts dominating. It is this behavior we see in figure 8: the computational time in the Monte Carlo method scales according to the first term, i.e., as L 2 .
Hence, we summarize: The time stepping procedure scales as L 3.99 whereas the Monte Carlo algorithm scales as L 1.98 , as shown in figure 8.

Limitations
A closer inspection of figures 4 to 7 shows that the match between the Monte Carlo and the time stepping procedures is good but not perfect. In this section we discuss the discrepancies between the two methods quantitatively. We show in figure 9 the non-wetting fractional flow for a 40 × 40 network using both time stepping and Monte Carlo with sub network size Λ ranging from 4 to 40. Notice that we also consider the sub-network size 40 which is equal to L. The calculations here are done in the constant Q ensemble with a capillary number Ca equal to 0.1 or 0.05. As we see, there is a systematic deviation between the time stepping and the Monte Carlo results that increases with increasing non-wetting saturation S nw . This deviation is highlighted in figure 10 where the difference between the time stepping and the Monte Carlo results for different Λ is shown. We note that the difference between the Monte Carlo and the time stepping decreases with increasing capillary number Ca. This is, however, to be expected, as for infinite Ca, any curve, Monte Carlo or time stepping, must fall on the diagonal of figure 9. In figure 11 we show the discrepancy between the pressure drop ∆P using time stepping and Monte Carlo for different sub lattice size Λ. The systematics seen in the fractional flow data, figures 9 and 10, where the difference grows with increasing non-wetting saturation is much less pronounced in this case.
In figure 12, we show histograms over the non-wetting saturation of the links. That is, we measure how much non-wetting fluid each link contains. When the overall non-wetting saturation S nw = 0.3, there is essentially no difference between the time stepping and the Monte Carlo result. However, for S nw = 0.8, there is a difference that depends on the sub lattice size Λ. This difference, measured as the area between the time stepping and the Monte Carlo histograms, is shown in figure 13 as a function of S nw . The picture seen here resembles that seen for the non-wetting fractional flow (figure 9): the difference grows with increasing S nw .
When the non-wetting saturation S nw is small, the non-wetting fluid will form bubbles or small clusters surrounded by the wetting fluid. As S nw is increased, these clusters grow in size until there is a percolation-type transition where the wetting fluid starts forming clusters surrounded by the non-wetting fluid. This scenario has been studied experimentally by Tallakstad et al. [23,24]. They argued that there is a length scale l * . Clusters that are larger than this length scale will move, whereas clusters that are smaller will be held in place by the capillary forces. The Monte Carlo algorithm calls for selecting a sub network which is then "lifted" out of the system, "folded" into a torus and then time stepped. The boundaries of the sub network will cut through clusters and mobilize these. This changes the cluster structure from that of the time stepping procedure. In order to investigate this we have studied the cluster structure in the model under Monte Carlo and time stepping. In order to do this, we identify the nonwetting clusters. To do this, two nodes are considered to be part of the same cluster if the link between them has a non-wetting saturation more than a threshold value, a clip-threshold c t . Here we use a clip threshold equal to c t = 0.9 [25]. In figure  14, we show typical cluster structures for two different non-wetting saturations obtained with Monte Carlo and with time stepping. For S nw = 0.7, the nonwetting clusters are still quite small and there is no discernable difference between the configurations obtained with time stepping and with Monte Carlo. However, for S nw = 0.8, there is one dominating cluster in the time stepping case whereas the clusters are more broken up in the Monte Carlo case.
We measure this qualitative difference in cluster structure for S nw = 0.8 by recording the cluster size distribution for the two types of updating, see figure  15. When following the time stepping procedure, we run the system for 500 pore volumes. During the last 125 pore volumes injected (1/4th of the total), we measure the cluster size distribution after passing each pore volume of fluids. When using Monte Carlo, we run the system for 400 Monte Carlo updates. We record the cluster size distribution for every of the last 100 updates. In both the time stepping and Monte Carlo runs, we average over 10 samples. The number of links belong to a cluster defines the size of that cluster. The total number of clusters is N total and the number of clusters of size k that we record is N k . We show P (k) = N k /N total in the figure. For S nw = 0.6 and 0.7, there is no discernable difference in the cluster structure between the Monte Carlo and the time stepping procedures. However, for S nw = 0.8, there are differences. For every k the number of clusters during the Monte Carlo updating procedure is larger than for the time stepping procedure, except for the largest clusters, the percolating cluster seen in figure  14. This supports the supposition that the Monte Carlo breaks up the large nonwetting clusters.
Clearly, for the Monte Carlo algorithm to be perfected, this tendency of chopping up large non-wetting clusters needs to be counteracted. Presumably, this is a problem that decreases with increasing system and sub lattice size as it is a boundary effect.

Conclusion
We have in this work presented a new Monte Carlo algorithm for immiscible twophase flow in porous media under steady-state conditions using network models. It is based on the Metropolis transition probability (23) which in turn is build upon the configuration probability (19) which we derive here. By steady-state conditions, we mean that the macroscopic parameters that describe the flow such as pressure difference, flow rate, fractional flow rate and saturation all have well defined means that stay constant. On the pore level, however, clusters flow, merge, break up, and so on. The flow may be anything but stationary. We described the algorithm in Section 3.2.2. Computationally, the Monte Carlo algorithm is very fast compared to time stepping. We find that the time stepping procedure when implemented on a square lattice demands a computing time that scales as the linear size of the lattice, L, to the fourth power, whereas the Monte Carlo method scales as the linear size to the second power, see Section 4.3. However, there is another term that contributes to the computing time in the Monte Carlo procedure which scales as L 4.88 . This term has a prefactor associated with it which is very small compared to the other term scaling as L 2 . For L up to about 230, this term is small compared to the first one.

Open Questions
There are open questions with respect to the Metropolis Monte Carlo approach that we present here. The most important step in the direction of constructing such an approach is to identify the configuration probability (19). The second most important step is to provide a way to generate trial configurations that obey the symmetry requirement (21). Section 3.2.2 is concerned with this.
We see three challenges that will need to be overcome before the Monte Carlo algorithm that we propose here is fully capable of replacing time stepping.
-The Monte Carlo algorithm needs to be generalized to irregular networks, e.g., those based on reconstructed porous media [1].
-The necessity to solve the Kirchhoff equations for the entire pore network once for every Monte Carlo update will slow down the algorithm when it is implemented for large systems. Ideally, one should find a way to circumvent this necessity.
-The Monte Carlo algorithm has a tendency to break up large non-wetting clusters as described in Section 4.4. This is a problem for large non-wetting saturations. It is most probably a boundary effect that comes from the way the sub networks are constructed. However, it needs to be overcome if the algorithm is to be useful for the entire range of saturations.
Overcoming these three challenges will allow network models to take advantage to the full of the ongoing revolution in pore space characterization.
We have in this article presented a first attempt at constructing a Markov Chain Monte Carlo algorithm based on the configurational probability (19). There is no reason not to believe that other ways of constructing such Monte Carlo algorithms might be possible that are both faster and do not pose the challenges listed above.