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

Abstract

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.

Keywords

1 Introduction

The characterization of porous media at the pore level is undergoing a revolution (Blunt et al. 2013). 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 scales—e.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 toward 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 toward this goal.

It is the aim of this paper to introduce a new algorithm that improves significantly the efficiency of network models (Joekar-Niasar and Hassanizadeh 2012). 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 is 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 bottleneck in this approach is the necessity to solve the Kirchhoff equations to determine the pressure field whose gradient drive the fluids in competitions with the capillary forces.

A different and at present popular computational approach, among several, is the lattice Boltzmann method (Ramstad et al. 2010, 2012). 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 matters, but not at the larger scale where the large-scale topology of the pore network is more important. Further methods that resolve the flow at the pore level are, e.g., smoothed particle hydrodynamics (Tartakovsky and Meakin 2005; Ovaysi and Piri 2010; Liu and Liu 2010) and density functional hydrodynamics (Armstrong et al. 2016). 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 scales 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 Sect. 2, we describe the network model we use to compare the Monte Carlo method with time stepping, see Aker et al. (1998), Knudsen et al. (2002). In the next Sect. 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 (Hansen et al. 2016; Savani et al. 2016). 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 Sect. 4 where we compare the Monte Carlo method with time stepping using the same network model described in Sect. 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 toward how these may be overcome. We end by Sect. 5 where we summarize the work and draw our conclusions.

2 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 Aker et al. (1998), Knudsen et al. (2002). For simplicity, we do not consider a reconstructed pore network based on a real porous medium (Tørå et al. 2011; Ramstad et al. 2009). Rather, we simply use a two-dimensional square network, with disorder in the pore radii, oriented at \(45^\circ \) with respect to the average flow direction as shown in Fig. 1. As described in Knudsen et al. (2002), 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\times 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=1\) mm. We neglect gravitational effects.

Geometry of the pore network we use. The shaded area constitutes a link between two nodes

Fluid flow through each link in the network is modeled using the Washburn equation (Washburn 1921), see Fig. 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_\mathrm{c}(x)\) where \(x \in [0,l]\) is the position of the interface. The capillary pressure is given by the Young–Laplace equation

where \(\gamma \) is the surface tension, \(\theta \) is the contact angle between the interface and the pore wall. We set \(\gamma \cos \theta =30\) dyn/cm. \(r_0\) is the average link radius. We assume that the link has a shape so that \(p_\mathrm{c}\) attains the given x dependence. It has been chosen so that \(p_\mathrm{c}(0)=p_c(l)=0\) and \(\max _x |p_\mathrm{c}(x)|=|p_\mathrm{c}(l/2)|\). The Washburn equation then becomes

where \(\mu _\mathrm{av}=s_{nw}\mu _{nw}+s_w\mu _w\) is the effective 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 \(\mu _{nw}= \mu _w=1\) poise.

This is one of the links in the network. The wetting and non-wetting fluids, colored by white and gray, respectively, are separated by interfaces. Each interface provides a capillary pressure \(p_\mathrm{c}(x)\) that points in the direction from the non-wetting toward the wetting fluid. Through the link, a flow q passes. We indicate the two node pressures \(p_1\) and \(p_2\) at the end of the link

A pressure difference \(\Delta 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. (2002). 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 corresponding matrix inversion problem by using the conjugate gradient algorithm (Batrouni and Hansen 1988).

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 Eq. (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 the 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. (2002).

3 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 (Press et al. 2007) after having presented a short review of the Metropolis version of Monte Carlo (Press et al. 2007; Metropolis et al. 1953).

3.1 Statistical Mechanics of Immiscible Two-Phase Flow

Sinha et al. (2013) 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\pi x/l]\) where \(l\ll L, a\) is an amplitude and \(r_0\) is the average radius of the tube. Suppose furthermore that the tube is filled with wetting fluid except for a bubble of length \(\Delta x_\mathrm{b}\) and a center position \(x_\mathrm{b}\). By using Eq. (1), one derives the net capillary force from the two interfaces that limit the bubble as,

where \(\pi r_0^2\dot{x}_\mathrm{b}=q\), and \(\Delta p=(l{/}L)\Delta P\) where \(\Delta P\) is the pressure difference across the tube.

Suppose there is a quantity \(f=f(x_\mathrm{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 \(\Delta P\) does not vary in time. The time average of f is then

where \(x_\mathrm{b}(t)\) is the time integration of the Washburn equation (2) and the time period \(T_\mathrm{b}=(2\pi \sigma )/\sqrt{\Delta p-\sigma ^2}\). We note, and this is the crucial observation, that we may change integration variable from time t to bubble position \(x_b\),

is the configuration probability. That is, the configuration of the tube is given by the position of \(x_\mathrm{b}\) of the bubble. Equation (8) gives the probability density to find the bubble at position \(x_\mathrm{b}\)—and hence in that configuration.

The Washburn equation (5) gives the motion of the bubble that is used in Eqs. (7) and (8). The Washburn equation assumes that we control the pressure drop \(\Delta P\). If we on the other hand control the flow q, the equation of motion becomes

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\in [1,N]\) are \(x_j\), and it has a width of \(\Delta x_j\). Since the system is one dimensional, all bubbles move with the same speed \(\dot{x}_j=\dot{x}_1\). The Washburn equation is then

Solving the equations of motion (12) gives \(x_j=x_j(t)\). We may invert \(x_1=x_1(t)\) to get \(t=t(x_1)\). Hence, we then have \(x_j(x_1)=x_j(t(x_1))\) for all j. Suppose now we have a function \(f=f(x_1,\ldots , x_N)\), analogous to the one introduced in Eq. (7). Its time average is

where \(q=\pi r_0^2 \dot{x}_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 denote that \({\mathbf {x}}=(x_1,x_1,x_2,\ldots ,x_N)\), where \(x_i\) is the position of the ith interface. 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 is highly nonlinear, but of the form \(\dot{x}_i=g_i({\mathbf {x}})\). Solving these equations gives \(x_j=x_j(t)\).

Again we consider a function \(f=f({\mathbf {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 \({\mathbf {x}}(t)={\mathbf {x}}(t(x_i))={\mathbf {x}}(x_i)\). 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 \(\dot{x_1}=Q/\pi r_0^2\), 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 \({\mathbf {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, Eq. (19) makes sense: The slower the flow, proportionally the more the system stays in—or close to—a given configuration (Savani et al. 2016).

Is the system ergodic? Equations (7), (14) and (16) answer this question positively. Time averages give, by construction, the same results as configurational averages.

3.2 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 (Krauth 2006; Landau and Binder 2015).

3.2.1 General Considerations

We have a set of configurations characterized by the variable \({\mathbf {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 \(\mathrm{d}{{\mathbf {x}}}\) of the configuration—is proportional to \(\varPi ({{\mathbf {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 \({\mathbf {x}}\) to state \({\mathbf {x}}'\) is constructed as

where \(\pi ({\mathbf {x}},{\mathbf {x}}')\) is the probability density to pick trial configuration \({\mathbf {x}}'\), given that the system is in configuration \({\mathbf {x}}\). It is crucial that \(\pi ({\mathbf {x}}', {\mathbf {x}})\) is symmetric,

3.2.2 The Implementation

The Metropolis Monte Carlo algorithm based on Eq. (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 \(\pi ({\mathbf {x}},{\mathbf {x}}')\), which must obey the symmetry (21). That is, if the system is in configuration \({\mathbf {x}}\), the probability to pick a trial configuration \({\mathbf {x}}'\) must be equal to the probability to pick as trial configuration \({\mathbf {x}}\) if the system is in configuration \({\mathbf {x}}'\).

Suppose the system is in configuration \({\mathbf {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 Fig. 3 and consists of \(L\times L\) links as described in Sect. 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

We show here a typical network of the kind we use for comparing the time stepping and Monte Carlo methods. The network is bi-periodic, and the flow is from the bottom toward the top. The dark red constitutes the non-wetting fluid and the gray constitutes the wetting fluid. When using the Monte Carlo method, a random sub-network is chosen as shown in the box, taken out of the network, integrated forward in time after having been made bi-periodic, and then re-entered into the network. This is the heart of the Metropolis Monte Carlo algorithm

We choose a randomly positioned sub-network as shown in Fig. 3. The network consists of \(\varLambda \times \varLambda \) 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 remain unchanged at this point.

By solving the Kirchhoff equations on the sub-network, we time step the configuration forwards in time while keeping the flow rate \(\varTheta \) 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 are 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 \({\mathbf {x}}'\).

Part of the probabilistic choice of the trial configuration that defines \(\pi ({\mathbf {x}},{\mathbf {x}}')\) rests on the choice of the sub-network: Its position is picked at random. Hence, if the system is in state \({\mathbf {x}}\) or in trial state \({\mathbf {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 \(\varTheta \). Hence, all sub-network configurations are equally probable, see Eq. (19). Hence, also this part of the choice of trial configuration is symmetric. The full probability \(\pi ({\mathbf {x}},{\mathbf {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 \(\varPi ({\mathbf {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 \({\mathbf {x}}\) has been generated, it is necessary to calculate the total flow rate \(Q=Q({\mathbf {x}}')\) in the network. We then decide to accept the trial configuration \({\mathbf {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.

4 Results

We now present numerical results of the Monte Carlo simulation considering the model described in Sect. 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 \(\Delta P\) is kept constant (CP ensemble). A network of \(40\times 40\) links (\(L=40\)) is considered for both Monte Carlo and time-stepping procedure. The sub-network size is \(20\times 20\) links (\(\varLambda =20\)). To identify whether the system has reached the steady state, we measured the quantities as a function of time steps in time stepping and as a function of sweeps in the case of Monte Carlo. We then identified the steady states when the averages of measured quantities (e.g., \(F_\text {nw}\) and \(\Delta P\) or Q) did not drift with time or with sweeps. We then take average over time (time stepping) or sweeps (Monte Carlo), which give us the time average and the ensemble average, respectively. We average 10 different networks, but with the same sequence of networks for both Monte Carlo and time stepping. First we present the results for CQ ensemble. Two capillary numbers, \(\text {Ca} = 0.1\) and 0.01, are used, and for each \(\text {Ca}\), simulations are performed for different values of non-wetting saturations in intervals of 0.05 from 0.05 to 0.95.

In Fig. 4, we plot \(F_\text {nw}\)—the non-wetting fractional flow—as a function of \(S_\text {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_\text {nw}=S_\text {nw}\), shown in the figure. Rather, \(F_\text {nw}\) is less than \(S_\text {nw}\) for low values of \(S_\text {nw}\) and higher than \(S_\text {nw}\) for higher value of \(S_\text {nw}\). It therefore crosses the \(F_\text {nw}=S_\text {nw}\) line at some point, which is not at \(S_\text {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 behavior is more prominent for the lower value of \(\text {Ca}\), as capillary forces play a more dominant role. The curves from the Monte Carlo and time-stepping calculations fall on top of each other for most of the lower to intermediate range of the saturation values, and we only see some difference at very high or low \(S_\text {nw}\). We will present a more quantitative comparison between the results of Monte Carlo and time stepping later in Sect. 4.4. The variation of total pressure drop \(\Delta P\) for the two capillary numbers as a function of \(S_\text {nw}\) is shown in Fig. 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_\text {nw}\). We only see differences at high values of \(S_\text {nw}\). \(\Delta P\) increases with \(S_\text {nw}\), reaching a maximum at some intermediate saturation and then decreases again. When \(S_\text {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 \(\Delta P\) after the maximum is due to the decrease in the number of interfaces blocking the fluids.

Values of pressure difference (\(\Delta P\)) for constant flow rate (Q) in the steady state as a function of non-wetting saturation (\(S_\text {nw}\)) for the capillary numbers \(\text {Ca}=0.1\) and 0.01 obtained via Monte Carlo (MC) simulations and time stepping (TS). The data are averaged over 10 samples

4.2 Constant \(\Delta P\) Ensemble

We now turn to the constant pressure ensemble. Here we keep \(\Delta P\) constant throughout the calculations. In this case, the Metropolis Monte Carlo algorithm, Eq. (23), becomes

Results for the simulations with constant \(\Delta P\) are shown in Figs. 6 and 7. Simulations are performed for two different values of \(\Delta P, 15\) and \(6.5\,\text {kPa}\). The steady-state values of \(F_\text {nw}\) show similar variation with \(S_\text {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_\text {nw}\). Here Q varies with the saturation, and the corresponding capillary numbers are plotted in Fig. 7 for Monte Carlo and time stepping. As discussed before, the number of interfaces first increases with the increase in saturation from zero, reaches a maximum value, and then decreases again as \(S_\text {nw}\) approaches 1. The pressure is constant here, so the total flow rate decreases with increasing capillary barriers at the interfaces and correspondingly \(\text {Ca}\) varies as in Fig. 7. Here again, good match between the results Monte Carlo and time stepping can be observed.

Non-wetting fractional flow (\(F_\text {nw}\)) as a function of non-wetting saturation in the steady state for constant \(\Delta P\) ensemble. Results are presented for Monte Carlo (MC) and time stepping (TS) for two different overall pressure drops \(\Delta P = 15\) and \(6.5\,\text {kPa}\). As Q varies with saturation for constant \(\Delta P, \text {Ca}\) is not constant here, which is demonstrated in the next Fig.7. The data are averaged over 10 samples

Capillary numbers, calculated from the total flow rates (Q), in the steady state as a function of the non-wetting saturation \(S_\text {nw}\) for constant \(\Delta P\) ensemble. Results are compared between Monte Carlo (MC) and time stepping (TS). The data are averaged over 10 samples

We show in Table 1 the percentage of rejections for the data shown in Fig. 7. The number of rejections is in all cases quite small. This can be understood as follows. Set \(Q({\mathbf {x}},\Delta P)=Q\) and \(Q({\mathbf {x}}',\Delta P)=Q+\delta \) where \(\delta \) may be positive or negative. Hence, the probability to accept the new configuration is

where we have assumed \(\delta \ll Q\). With a small \(\delta \), the probability to reject the trial configuration is small. This is reflected in Table 1.

Table 1

Percentage of rejected configurations in the constant \(\Delta P\) ensemble

\(\Delta p\)

\(S_{nw}\)

Rejections (%)

15 kPa

0.3

2.1

0.5

2.3

0.7

1.5

6.5 kPa

0.3

8.8

0.5

11.6

0.7

4.2

4.3 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_\text {MC}\) for the Monte Carlo method and \(T_\text {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\times 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^\beta \), where \(\beta =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, \(\beta \) is much smaller than four.

The number of time steps needed to push one pore volume through the network is \(n_\text {t}\). We expect it to depend on L as \(n_\text {t}=aL^\tau \), 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 length of the network, L, making \(\tau =1\). In practice, as we shall see, it is slightly larger.

For each time step, the conjugate gradient demands \(t_\text {cg}=b L^\beta \) operations where b is another prefactor. The total computational time (\(T_\text {TS}\)) per pore volume is then

where \(\alpha _\text {TS}=\tau +\beta \). Based on the theoretical considerations above, setting \(\beta =4\) and \(\tau =1\), we have \(T_\text {TS}\sim L^5\). The actual computational time measured using the clock() function in C is plotted in Fig. 8 for \(\text {Ca}=0.1\) and \(S_\text {nw}=0.4\). We find that \(T_\text {TS}\) scales with L with an exponent \(\alpha _\text {TS}=3.99\pm 0.03\), which is much smaller than 5. Measuring \(n_\text {t}\) and \(t_\text {cg}\) independently gives \(\tau =1.11\pm 0.03\) and \(\beta =2.88\pm 0.02\), see the inset in Fig. 8.

Plot of total computational time, \(T_\text {MC,TS}\) (in seconds) used by Monte Carlo (MC) and time stepping (TS) for different system sizes (L). 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 \(\varLambda =L=20\), the timing of the two methods is equal. We use the CQ ensemble with \(\text {Ca}=0.1\) and \(S_\text {nw}=0.4\). Six different system sizes, \(L=20, 40, 60, 80\) and 120, are considered. From the slopes, the exponents \(\alpha \) for time stepping and Monte Carlo are found. For Monte Carlo, we find \(\alpha _\text {MC}= 1.98\pm 0.01\), which is close to the theoretically expected value \(T_\text {MC}\sim L^2\) (see text). However, for time stepping, we find \(\alpha _\text {TS} = 3.99\pm 0.03\), which is much smaller than theoretical expectation—\(T_\text {TS}\sim L^5\). In the inset, we plot the average time, \(t_\text {cg}\), taken by the conjugate gradient solver to solve one entire pressure field. We find \(t_\text {cg}\sim L^{2.88\pm 0.02}\). The number of time steps per pore volume, \(n_\text {t}\), scales as \(n_\text {t}\sim L^{1.11\pm 0.03}\). Combining these two results, we find that the computational time for the time-stepping procedure to scales as \(T_\text {TS}\sim L^{3.99}\)

For the Monte Carlo algorithm, each sweep ideally contains \((L/\varLambda )^2\) individual Monte Carlo updates. Each Monte Carlo update consists of time stepping a sub-lattice of size \(\varLambda \times \varLambda \). Hence, the cost of a Monte Carlo update is \(ab\varLambda ^{\alpha _\text {TS}}\) when using Eq. (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 \(b L^\beta \). The time per Monte Carlo sweep is then

where \(\alpha _\text {TS}=3.99\) and \(\beta =2.88\). The factor “4” signifies that we time step the sub-lattice for four pore volumes. By setting \(a\approx 10\) and \(\varLambda =20\), the first term will dominate compared to the second term on the right-hand side of this equation if \(4a\varLambda ^{\alpha _\text {TS}} \approx 6.4\times 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 Fig. 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 Fig. 8.

4.4 Limitations

A closer inspection of Figs. 4, 5, 6 and 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 Fig. 9 the non-wetting fractional flow for a \(40\times 40\) network using both time stepping and Monte Carlo with sub-network size \(\varLambda \) 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_\text {nw}\). This deviation is highlighted in Fig. 10 where the difference between the time stepping and the Monte Carlo results for different \(\varLambda \) 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 \(\text {Ca}\), any curve, Monte Carlo or time stepping, must fall on the diagonal of Fig. 9.

Non-wetting fractional flow \(F_\text {nw}\) as a function of non-wetting saturation \(S_\text {nw}\) for time stepping compared to Monte Carlo for different sub-network sizes (\(\varLambda \)) in the constant Q ensemble. The size of the network, L, is 40 for both Monte Carlo and time stepping

Difference of the non-wetting fractional flow (\(\Delta F_\text {nw}\)) between time stepping and Monte Carlo for different values of \(\varLambda \) is plotted as a function of \(S_\text {nw}\). \(\Delta F_\text {nw}\) fluctuates around zero for \(\varLambda = L\), and a systematic increase is observed with the decrease in \(\varLambda \) for the whole range of \(S_\text {nw}\)

In Fig. 11, we show the discrepancy between the pressure drop \(\Delta P\) using time stepping and Monte Carlo for different sub-lattice size \(\varLambda \). The systematics seen in the fractional flow data, Figs. 9 and 10, where the difference grows with increasing non-wetting saturation is much less pronounced in this case.

Pressure difference \(\Delta P\) as a function of non-wetting saturation \(S_\text {nw}\) for time stepping compared with Monte Carlo for different sub-network sizes (\(\varLambda \)) in the CQ ensemble. The size of the network, L, is 40 for both Monte Carlo and time stepping

In Fig. 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_\text {nw}=0.3\), there is essentially no difference between the time stepping and the Monte Carlo result. However, for \(S_\text {nw}=0.8\), there is a difference that depends on the sub-lattice size \(\varLambda \). This difference, measured as the area between the time stepping and the Monte Carlo histograms, is shown in Fig. 13 as a function of \(S_\text {nw}\). The picture seen here resembles that seen for the non-wetting fractional flow (Fig. 9): the difference grows with increasing \(S_\text {nw}\).

Comparison of cumulative distribution \(P(s_i)\) of link saturation \(s_i\) for time stepping and for Monte Carlo with different sub-system sizes. For \(S_\text {nw}=0.3, P(s)\) for Monte Carlo match with time stepping for all the sub-system sizes, whereas for \(S_\text {nw}=0.8\), a systematic difference in P(s) is observed for \(\varLambda <L\)

Area between the \(P(s_i)\) curves (Fig. 12) for time stepping and for Monte Carlo with different sub-system sizes as a function of \(S_\text {nw}\)

When the non-wetting saturation \(S_\text {nw}\) is small, the non-wetting fluid will form bubbles or small clusters surrounded by the wetting fluid. As \(S_\text {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. (2009a, b). 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 non-wetting 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_\mathrm{t}\). Here, we use a clip threshold equal to \(c_\mathrm{t} = 0.9\) (Ramstad and Hansen 2009). In Fig. 14, we show typical cluster structures for two different non-wetting saturations obtained with Monte Carlo and with time stepping. For \(S_\text {nw}=0.7\), the non-wetting 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_\text {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.

Typical non-wetting clusters for Monte Carlo (MC) and time stepping (TS) at \(\text {Ca}=0.05\). The network is of \(40\times 40\) links, and the sub-network size for Monte Carlo is \(20\times 20\) links. Each cluster is marked with different colors so that the structure is readily visible

We measure this qualitative difference in cluster structure for \(S_\text {nw}=0.8\) by recording the cluster size distribution for the two types of updating, see Fig. 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_\mathrm{total}\) and the number of clusters of size k that we record is \(N_k\). We show \(P(k)=N_k/N_\mathrm{total}\) in the figure. For \(S_\text {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_\text {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 Fig. 14. This supports the supposition that the Monte Carlo breaks up the large non-wetting clusters.

Cluster size distribution \(P(k)=N_k/N_\mathrm{total}\) versus cluster size k for time stepping and Monte Carlo. The blue circles signify the Monte Carlo data and the red circles the time-stepping data. The red and blue curves with triangles pointing upwards or downwards signify the Monte Carlo and time-stepping data after logarithmic binning. Here \(L=40\) and \(\varLambda =20\). The data are averaged over 10 samples

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.

5 Conclusion

We have in this work presented a new Monte Carlo algorithm for immiscible two-phase 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 Sect. 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 Sect. 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.

5.1 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 (Blunt et al. 2013).

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 Sect. 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.

Notes

Acknowledgements

IS, AH and SK thank VISTA, a collaboration between Statoil and the Norwegian Academy of Science and Letters for financial support. SS and AH thank the Beijing Computational Science Research Center and its director, Professor Hai-Qing Lin, for financial support and for providing an excellent atmosphere for doing science. We thank Eirik Grude Flekkøy, Knut Jørgen Måløy, Miguel Rubi and Marios Valavanides for many interesting discussions.

Copyright information

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.