Pressure and flow statistics of Darcy flow from simulated annealing

TL;DR

This paper demonstrates that simulated annealing can effectively compute pressure and flow statistics in Darcy flow through random media, matching traditional methods and offering advantages for complex permeability distributions.

Contribution

It introduces a novel application of simulated annealing for Darcy flow statistics, highlighting its flexibility and computational benefits over conventional methods.

Findings

Simulated annealing results agree with finite-volume calculations.

The method handles arbitrary permeability distributions.

Advantages include no matrix inversion and suitability for multiphase flow.

Abstract

The pressure and flow statistics of Darcy flow through a random permeable medium are expressed in a form suitable for evaluation by the method of simulated annealing. There are several attractive aspects to using simulated annealing: (i) any probability distribution can be used for the permeability, (ii) there is no need to invert the transmissibility matrix which, while not a factor for single-phase flow, offers distinct advantages for the case of multiphase flow, and (iii) the action used for simulated annealing is eminently suitable for coarse graining by integrating over the short-wavelength degrees of freedom. In this paper, we show that the pressure and flow statistics obtained by simulated annealing are in excellent agreement with the more conventional finite-volume calculations.

Figures

Tables

Table 1. Table 1: Values used for the variance σ 2 superscript 𝜎 2 \sigma^{2} of the log-permeability.

No.	Log-permeability variance ${\sigma }^2$	Color
1	$0.125$	gray
2	$0.25$	red
3	$0.5$	blue
4	$1.0$	magenta
5	$1.75$	brown
6	$2.5$	orange

Equations

$$\mathbf{q}=-K\nabla \mathbf{p}.$$

$$\nabla \cdot \mathbf{q}=0.$$

$${d\over dx}\biggl{(}K{dp\over dx}\biggr{)}=0\,,$$

$$K\frac{dp}{dx}=q_0,$$

$$\nabla \cdot (K\nabla p)=0,$$

$$K\nabla p=\nabla \times \mathbf{f}+\nabla g,$$

$$K(\mathbf{x})=e^{L(\mathbf{x})},$$

$$\displaystyle Z=\int\prod_{i}dp_{i}\int\prod_{i}{dL_{i}}\exp\biggl{[}-\sum_{ij}L_{i}({\mathsf{C}}_{L}^{-1})_{ij}L_{j}\biggr{]}$$

$$\displaystyle\quad\times e^{-\sum_{i}L_{i}}\delta\biggl{(}{p_{i}-p_{i-1}\over\Delta x}+q_{0}e^{-L_{i}}\biggr{)}\,.$$

$$\exp\biggl{[}-\sum_{ij}L_{i}({\mathsf{C}}_{L}^{-1})_{ij}L_{j}\biggr{]}$$

$$Q(\{p_i\})=\frac{e^{-S(\{p_i\})}}{Z},$$

$$\displaystyle S(\{p_{i}\})=\sum_{i}\log\biggl{(}{p_{i-1}-p_{i}\over q_{0}\Delta x}\biggr{)}$$

$$\displaystyle+\sum_{ij}\log\biggl{(}{p_{i-1}-p_{i}\over q_{0}\Delta x}\biggr{)}({\mathsf{C}}_{L}^{-1})_{ij}\log\biggl{(}{p_{j-1}-p_{j}\over q_{0}\Delta x}\biggr{)}\,.$$

$$\displaystyle Z_{\mathrm{2D}}=\int\prod_{ij}dp_{ij}\int\prod_{kl}{dL_{kl}}\exp\bigg{(}-\sum_{ij}L_{ij}\bigg{)}$$

$$\displaystyle\quad\times\exp\bigg{[}-\sum_{ij,kl}L_{ij}({\mathsf{C}}_{L}^{-1})_{ij,kl}L_{kl}\biggr{]}$$

$$\times \delta \left\{\frac{\partial }{\partial x}\left(e^{L_{ij}}\frac{\partial p_{ij}}{\partial x}\right)+\frac{\partial }{\partial y}\left(e^{L_{ij}}\frac{\partial p_{ij}}{\partial y}\right)\right\}.$$

$$\frac{\partial }{\partial x}\left(e^{L_{ij}}\frac{\partial p_{ij}}{\partial x}\right)$$

$$=\frac{1}{\Delta x^2}\left[e^{L_{i,j}}{p}_{i,j}-(e^{L_{i,j}}+e^{L_{i-1,j}})p_{i-1,j}+e^{L_{i-1,j}}{p}_{i-2,j}\right]$$

$$\frac{\partial }{\partial y}\left(e^{L_{ij}}\frac{\partial p_{ij}}{\partial y}\right)$$

$$=\frac{1}{\Delta y^2}\left[e^{L_{i,j}}{p}_{i,j}-(e^{L_{i,j}}+e^{L_{i-1,j}})p_{i-1,j}+e^{L_{i-1,j}}{p}_{i-2,j}\right]$$

$$\displaystyle Z_{\mathrm{2D}}=\int\prod_{ij}dp_{ij}\int\prod_{kl}{dL_{kl}}\exp\bigg{(}$$

$$\displaystyle\quad\times\lim_{t\to 0}\bigg{\{}\exp\biggl{[}{1\over t}\biggl{(}$$

$$+$$

$$\sum_{ij}\biggl{\{}\frac{\partial}{\partial x}\left(e^{L_{ij}}\frac{\partial p_{ij}}{\partial x}\right)+\frac{\partial}{\partial y}\left(e^{L_{ij}}\frac{\partial p_{ij}}{\partial y}\right)\biggr{\}}.$$

$$S=\frac{1}{2}{\int }_{V}K(\mathbf{x}){\left(\nabla p(\mathbf{x})\right)}^{2}dV,$$

$$S[p(\mathbf{x})+\delta p(\mathbf{x})]-S[p(\mathbf{x})]=0.$$

$$S[p+\delta p]-S[p]=\int K\left(\nabla p\right)\left(\nabla \delta p\right)dV$$

$$=-\int \nabla \cdot \left(K\nabla p\right)\delta pdV-{\int }_{\partial V}K\delta p\nabla p\cdot d\mathbf{S}.$$

$$\nabla \cdot (K\nabla p)=0,$$

$$N_x\times N_y\times N_z=50\times 70\times 50$$

$$X\times Y\times Z=40\mathrm{m}\times 85\mathrm{m}\times 25\mathrm{m}.$$

$${\lambda }_x=8\mathrm{m};\text{leavevmode}{\lambda }_y=8\mathrm{m};\text{leavevmode}{\lambda }_z=5\mathrm{m}.$$

$$p_i^{\mathrm{new}}=2p_i^{\ast }-p_i^{\mathrm{old}}$$

$$T^{(k)}={\alpha }^k{T}_i,$$

$$\mathcal{O}(2N_x{N}_y{N}_z\log (N_x{N}_y{N}_z))$$

Full text

Darcy flow, random permeability, simulated annealing, pressure statistics, flow statistics

I Introduction

The study of flow in porous media has a wide range of applications, including hydrology (e.g. bear ; mualem ), geothermal engineering (e.g.abdelazizac ), materials science (e.g. bear ), and the medical sciences (e.g. khaled ). Another application, and the focus of the work reported here, is the flow of oil in a reservoir.

There are two basic ways to conceptualize the flow through a porous material. One approach is to solve the Navier-Stokes equations. On domains in the millimeter to centimeter scale, the fluid configurations can be imaged with X-ray microcomputed tomography liu17 . The flow of oil through a rock, by contrast, calls for flow descriptions on the kilometer scale. Solving the Navier–Stokes equations on the later scale is not feasible because of limited information about the rock permeability and the matrix form in which to cast the problem. But, even if all of this information was available, the computational requirements would be prohibitive.

On a more coarse grained scale, such as the slow flow of a viscous fluid, the Navier-Stokes equations can be reduced to Darcy’s law, a relation between the effective permeability $K$, the average velocity (“flow”) $\bf{q}$ of the fluid, and the pressure gradient $\nabla\bf{p}$:

[TABLE]

Here, $K({\bf x})=k({\bf x})/{\mu}$, where $k({\bf x})$ is the effective permeability of the medium and $\mu$ is the viscosity of the fluid. The effective permeability describes the medium on a “mesoscopic” scale, large compared to the pore scale, but small on the scale of the macroscopic medium. Although proposed as an empirical relation by Darcy in the 1850s darcy , Darcy’s law can be recovered from the Navier-Stokes equations whitaker . In this paper, we focus on single-phase, incompressible flow:

[TABLE]

In reservoir engineering, Darcy’s law is solved numerically using the finite-volume method aziz ; schafer , while the finite-element method is common in hydrology.

In previous work westbroek , we formulated the solution to one-dimensional Darcy flow as a path integral over pressure. In discrete form, the path integral is a tool to simulate Darcy pressure paths $\{p_{i}\}$ on a spatial lattice using Markov chain Monte Carlo methods according to the weighting $e^{-S}$, where the “action” $S=S[\{p_{i}\}]$ contains the solution to Darcy’s law.

In one-dimension, the path integral hinges on an analytical solution obtained by combining (1) and (2),

[TABLE]

from which we immediately obtain

[TABLE]

where $q_{0}$ is a constant. In higher dimensions, the equation corresponding to (3),

[TABLE]

has the general solution

[TABLE]

where ${\bf f}$ is a twice differentiable vector function and $\nabla^{2}g=0$. This solution is not as amenable to the path integral formulation as (4), so we turn to an alternative expression for the action that utilizes simulated annealing kirkpatrick ; press .

Although solving the problem using simulated annealing is more computationally intensive than the finite-volume method, there are several advantages over various other techniques. Simulated annealing allows any probability distribution to be used for the permeability, as does the finite-volume method. However, there is no need to invert the transmissibility matrix during simulated annealing, which, while not an issue single-phase flow, is advantageous for multiphase flow. On a more formal level, the action used for simulated annealing is eminently suitable for integration over the short-range degrees of freedom to derive coarse-grained permeability coefficients hristopulos ; attinger ; eberhard ; hanasoge ; king2 . The porous medium can also be characterized as a statistically homogeneous continuum with local fluctuations in physical parameters. The resulting path integral expression can be averaged over parameter fluctuations to obtain large-distance parameters that describe the flow tanksley .

Our paper is organized as follows. Path integral formulations of Darcy flow through random porous media are derived in Sec. II. We explain why the procedure we have developed previous for the numerical evaluation of path integrals in one dimension cannot be extended to higher dimensions. Section III discusses computation methods, including the simulation of the permeability fields and the method of simulated annealing. Results are presented in Sec. IV. We have calculated empirical probability densities for the pressure, Cartesian flow components, the total flow along the $y$-direction, and examined the effect of the variance of the log-permeability. Conclusions and a discussion are provided in Sec. V, including an assessment of the viability of simulated annealing and the extension of the path integral approach to multi-phase flow.

II Theory

II.1 Effective permeability

The effective permeability is modelled as a stochastic process. Various such models exist, with an overview is given in renard . We have opted to model $K({\bf x})$ as a lognormal process,

[TABLE]

where $L({\bf x})$ is a Gaussian process with mean zero, characterized by its variance $\sigma^{2}$ and correlation length $\xi$. More information on this process can be found in long1D . This conventional choice law has the advantage of a strictly positive permeability. We emphasize, however, that our method applies to any choice of permeability distribution.

II.2 Path integral in higher dimensions

In an earlier paper westbroek , we developed a path integral for Darcy flow in one dimension, The integral is over all discrete pressure trajectories that (subject to the boundary conditions) follow Darcy’s law, which is enforced by a delta-functional:

[TABLE]

The term

[TABLE]

encodes the correlated Gaussian probability distribution of the log-permeability, where $C_{L}$ denotes the log-permeability covariance matrix. The factor $e^{-\sum_{i}L_{i}}$ is the Jacobian associated with integrating over the $L_{i}$ rather than over the $e^{L_{i}}$. Upon integration over the $L_{i}$ we obtain a path integral

[TABLE]

with the discrete “action” westbroek

[TABLE]

Discrete pressure paths are generated according to the probability density (10).

An analogous path integral to (II.2) in two dimensions is obtained with the standard procedure for classical statistical dynamics phythian77 ; peliti78 ; jouvet79 ; jensen81 :

[TABLE]

The delta-function enforces the discrete form of Darcy’s law (5), where we have used the notation

[TABLE]

The next step is to represent the delta-function as the limit of an exponential, so that the exponentials in (12) can be combined into a single exponential whose argument is the “action”. The usual procedure phythian77 ; peliti78 ; jouvet79 ; jensen81 is to apply a functional Fourier transform, which yields a complex action. This is appropriate for formal studies involving perturbation expansion, where the complex action yields real results despite the complex nature of intermediate calculations. But the Markov chain Monte Carlo method relies on real variables from the outset, so we represent the delta functional as the limit of a Gaussian probability density:

[TABLE]

This expression is readily generalized to three dimensions.

Averages of pressure and correlation functions can be calculated from (II.2) by first generating permeability fields, then setting $t$ to some value, and finally using the Metropolis–Hasting (MH) algorithm to minimize the discrete “action”:

[TABLE]

Successively smaller values of $t$ are chosen until there is convergence of the pressure distributions. This procedure, which requires separate calculations for each value of $t$, is not especially efficient. In the next section, we develop a more elegant approach based on simulated annealing.

II.3 Path integral for simulated annealing

To account for any number of dimensions, let us write the action in its continuum form:

[TABLE]

where the integral is carried out over the entire volume under consideration. Simulated annealing aims to find the pressure $p({\bf x})$ that minimizes the action (16) for fixed permeability $K({\bf x})$.

We show that the minimized pressure follows Darcy’s law. The objective is to extremize the action (16) with respect to the pressure. We vary the action with respect to $p$ by adding an infinitesimal pressure $\delta p$ and imposing the condition

[TABLE]

Any boundary conditions are unchanged, so $\delta p({\bf x})=0$ at the boundaries of the volume $V$. If there exists a $p^{\ast}$ such that the stationarity condition (17) holds, the action $S$ is stationary at $p^{\ast}$. Retaining $\delta p$ only to first order and performing an integration by parts, we obtain

[TABLE]

The boundary term vanishes because of Dirichlet boundary conditions fix the pressure across the entry and exit surfaces, and the absence of pressure fluctuations along surfaces perpendicular to the flow direction (Sec. III). Because $K$ is always nonzero, the condition (17) translates into

[TABLE]

which is Darcy’s law (5) for incompressible flow. Thus, if $p^{\ast}$ can be found such that the stationarity condition is met at fixed $K$, then $p^{\ast}$ follows Darcy’s law. The simulated annealing algorithm can be applied to the action (16) to solve for the pressure. The method is inspired by the process of annealing, which is a treatment whereby a solid is slowly cooled until its structure is eventually frozen at its minimum free energy configuration du .

III Computational methods

We will follow the convention for units used by the hydrology community. Darcy’s law (1) gives the relation between the effective permeability $K\leavevmode\nobreak\ \mathrm{[L^{2}\,T^{-1}]}$, the flow $\bf{q}\leavevmode\nobreak\ \mathrm{[L\,T^{-1}]}$ and the pressure (also known as the hydraulic head) $p\leavevmode\nobreak\ \mathrm{[L]}$. The total flow in a given direction is denoted by $Q_{i}\leavevmode\nobreak\ \mathrm{[L^{3}\,T^{-1}]}$ $(i=x,y,z)$. Our three-dimensional, rectangular prismatic porous medium is simulated on a grid of

[TABLE]

lattice elements, representing a domain of size

[TABLE]

The correlation lengths of the permeability field are

[TABLE]

We have used Dirichlet boundary conditions along the $y$-direction, making that the main flow direction. No-flow boundary conditions were imposed along the other boundaries. The values chosen for the log-permeability variance $\sigma^{2}$, the defining feature of the six parameter sets, are given in Table 1. The values in this table were those used by Nowak et al. nowak_main , which enables a qualitative comparison to be made between our two approaches.

III.1 Simulation of the permeability fields

To simulate the permeability field, the first step is to generate the three-dimensional log-permeability field. Once the Gaussian field $L({\bf x})$ is available, definition (7) can be invoked to calculate the permeability field $K({\bf x})$. To generate $L({\bf x})$, we have made use of the circulant embedding technique gneiting , where the correlation matrix ${\sf C}$ of the desired field is embedded into a matrix ${\sf M}$ that has a circulant or block circulant structure. Products of the square root ${\sf M}^{1/2}$ with white noise random vectors are realizations of the desired random field nowak ; dietrich . This method relies on the fast Fourier transform (FFT). For a $d$-dimensional rectangular mesh containing $N_{x}^{d}$ points, the computational requirements are those of an FFT of a vector of dimension $2N_{x}^{d}$ per realization dietrich .

III.2 Simulated annealing

The simulated annealing algorithm is based on the MH algorithm, a step-by-step explanation of which can be found in westbroek2 . In the present case, simulated annealing seeks to minimize the action (16). Clearly, the minimum attainable value is zero. The algorithm consists of the following steps.

Initialize a random pressure that is consistent with the boundary conditions. 2. 2.

Execute the MH algorithm some $M\gg 1$ times. The MH algorithm lowers the value of the action $S$, but also accepts some modifications to the pressure that increase the action. It explores the entire “state space” (the set of values of $S$ as a function of $p({\bf x})$) and does not get stuck in a local minimum of the state space. 3. 3.

After every $N_{s}\gg 1$ steps, check the value of $S$. When the value of $S$ starts to fluctuate around a constant value, go to step 4. 4. 4.

Adapt the MH accept/reject criterion to “accept the change in the action with probability $\mathrm{min}(1,e^{-\delta S/T})$ for some constant $0<T<1$”. This is a “cooling step” dimensionAction . The state space is explored in smaller steps than was the case for the standard MH algorithm, while maintaining a constant acceptance rate. The lower the value of $T$, the smaller the steps. In our context, $T$ does not have the interpretation of a temperature, but its effect remains that of slowing down the state space exploration. 5. 5.

Repeat steps 3-4 until the action attains a critically low value $\epsilon_{1}$, say, $\epsilon_{1}=0.1$. 6. 6.

Employ a modification of the MH algorithm known as the “greedy algorithm”, which accepts only changes to the pressure that lower the action, until $S$ dives below a second critical value $\epsilon_{2}$, say $\epsilon_{2}=10^{-2}$.

To expedite the simulated annealing algorithm, we made use of a technique known as over-relaxation (OR). The idea behind over-relaxation is to choose trial changes that cause significant changes to the pressure field, but only small changes to the action creutz2 ; brown . Such strategic updates combine a thorough exploration of the phase space with a high probability of acceptance. Because the action (16) is quadratic, it is possible to calculate an update that leaves the action unchanged, rendering the Metropolis accept/reject step unnecessary. The update

[TABLE]

lies “on the other side” of the minimum of the action: $p_{i}^{\ast}$ is the value of $p_{i}$ that minimizes $S$, with all other parameters kept fixed. Since the over-relaxation procedure is deterministic, we alternate between OR and regular MH steps to avoid any dependence of the pressure field on its random starting configuration. Here, we exchanged three in four Metropolis sweeps for OR sweeps.

We have used an exponential cooling scheme

[TABLE]

where $k$ indicates the cooling step. The MH algorithm was executed $M=2,000$ times for all parameters. For $\sigma^{2}\leq 1$, repeating the cooling algorithm $N_{s}=3,000$ times was found to be a good choice. For $\sigma^{2}>1$, it was necessary to set $N_{s}=6,000$. To be able to calculate empirical probability density functions for the pressure and flow, we have worked with $N=10,000$ realizations for each parameter set.

In comparing the computational efforts involved in running the FVM method and simulated annealing, we note that both require a permeability field as input. The computational cost associated with the FFT is

[TABLE]

floating point operations (“flops”). The key calculation in the FVM is a sparse matrix inversion. The sparse matrix solver UMFPACK UMFPACK can solve such an equation in

[TABLE]

flops. Contrary to the FVM, the simulated annealing requires $\mathcal{O}((N_{x}N_{y}N_{z})^{2})$ flops to calculate a pressure realization. One factor $N_{x}N_{y}N_{z}$ arises from the number of lattice sites. The number of required intermediate updates $N_{\mathrm{sep}}$ introduces a further factor $N_{x}N_{y}N_{z}$. However, there are techniques whose implementation is likely to decrease the run time considerably, such as the multigrid Monte Carlo (MGMC) method goodman ; janke ; liu and directed sampling duncan1 ; duncan2 ; leli .

IV Results

We have calculated empirical probability densities for the pressure $p$, the flow components $q_{y}$ and $q_{x}$, and the total flow $Q_{y}$ and looked at the effect of the variance of the log-permeability $\sigma^{2}$. Based on information about boundedness, we have made parametric fits using the choices made in nowak_main for guidance.

All quantities were normalized for straightforward interpretation. The Dirichlet boundary conditions were chosen to yield a pressure difference

[TABLE]

The flow components were normalized as

[TABLE]

where $K_{e}$ is the theoretical expectation value of the permeability, $I_{0}=\Delta p/L_{y}$ and $A=L_{x}\times L_{z}$. Computationally, the normalization was achieved by setting $K_{e}=1/(I_{0}A)$.

Due to correlations in the permeability, the distributions of pressure and flow are in general not expected to be Gaussian, especially when the variance of the log-permeability is high. The pressure in the main direction takes values in the interval $[0,1]$, due to the Dirichlet boundary conditions. Given this constraint and the choice of stochastic model for the permeability, the log-normal distribution is an obvious contender for parametric fits to the pressure. The probability density function is given by:

[TABLE]

In order to visualize the dependence of the pressure on its position along the main axis, we have made empirical probability density plots at two positions: $(0.5X,0.5Y,0.5Z)$, which is the center of the domain (Fig. 1) and $(0.5X,0.8Y,0.5Z)$ (Fig. 2). From the two sets of figures it is apparent that the lognormal distribution is most evident near the boundary. Towards the center of the domain, the histograms bear more resemblance to the normal distribution ababou , as the generalized Central Limit Theorem predicts bouchaud . For a more extensive explanation of this theorem in the context of Darcy flow, see long1D . These boundary effects increase with the log-permeability variance, as can also be observed for the flow.

Like the pressure, the flow along the main axis is subject to a non-negativity constraint, which enforces a one-sided bound. The flow could only be negative in the unlikely event of flow reversal due to locally very high permeability. The observed values for the cases considered in this work were non-negative. For the flow, as for the pressure, the shape of the probability density depends on the vicinity to a restricting boundary. We have evaluated the flow in the main direction at the center of the domain. The results are shown in Fig. 3.

One can see that for high values of the log-permeability variance, the flow statistics are most clearly log-normal. This is because the boundary effects are more strongly felt for high values of $\sigma^{2}$, a pattern that can also be observed by comparing Figures 1 and 2. The total flow in the main direction, defined as the average over a cross-section perpendicular to the $y$-axis, is conserved along the $y$-axis. The results can be seen in Fig. 4. When comparing Figs. 3 and 4, an obvious difference is that the flow statistics of the total flow approximate the Gaussian distribution more closely. The log-normal distribution tends to the normal distribution in the limit $(\sigma^{\prime}/\mu^{\prime})^{2}\to 0$. The Gaussian appearance is a result of the averaging over a cross-section that defines the total flow.

For the transverse flow $q_{x}^{*}$, shown in Fig. 5, we fitted an exponential power distribution. This choice reflects the expectation that the transverse flow is symmetric about zero. The probability density function for the exponential power law is:

[TABLE]

The Gaussian distribution corresponds to the case $k=2$. All parametric fits for the parameters $\{\mu^{\prime},\mu^{\prime\prime},\sigma^{\prime},\sigma^{\prime\prime},k\}$ were made using the in-built routine FindFit of Mathematica mathematica .

A striking feature of the transverse flow statistics are the long and heavy tails. The exponential power distribution was chosen because it reflects this feature. This choice was also made in nowak_main ). The tails are heaviest for high values of the log-permeability variance. In the limit of high variance, the permeability can take a very wide range of values. Thus, it will often occur that the flow either continues along its main axis, or is diverted. This behavior is reflected in the statistics by the tails of the distribution.

V Conclusion and Discussion

The rock heterogeneities exert a significant influence on the flow, from the pore scale up to the kilometer scale. Calculation of the Darcy pressure statistics depends on an explicit description of the permeability $K$ at the mesoscopic scale. For this work, we have chosen a lognormal distribution. Simulated annealing can be used to calculate the pressure and flow statistics for any type of realization of the permeability $K$. An alternative to assuming the lognormal distribution could be the use of multiple-point statistics, a method that directly infers the necessary multivariate distributions from training images huysmans , or copulas, which describe the stochastic structure without reference to the corresponding marginal distributions bardossy .

We have shown that our action $S$ (16) can be used to apply simulated annealing to calculate Darcy pressure and flow statistics. We have outlined our computational methods in such a way as to make them easily reproducible. Our model was a three-dimensional, bounded domain, with Dirichlet boundary conditions at two ends and no-flow boundary conditions at the remaining four. Our results for the pressure, calculated at two different points in the domain, as well as those for the local and total flow in the main and in a transverse direction al behaved qualitatively as expected. Parametric log-normal and exponential power-law fits were made using Mathematica, all of which passed one-sided Kolmogorov-Smirnov tests at the $95\%$ confidence level. At the moment, simulated annealing is not computationally competitive with the finite-volume method. Its runtime may be improved through the use of the multigrid method, however. Most promisingly, it may be possible to apply the renormalization group as an upscaling technique. Such an application would provide a different take on the problem and may enable the user to run coarse but fast simulations to capture the main characteristics of the pressure and flow statistics.

A major challenge is the extension of the present approach to multiphase flow. For two-phase flow a generalized form of Darcy’s law is used:

[TABLE]

where the subscript $i$ represents the fluid phase (oil or water), $k_{r,i}$ is the relative permeability, and $S_{i}$ is the pore volume fraction of the fluid phase $i$. The two volume fractions must sum to one. The total velocity is given by the sum of the individual phase velocities:

[TABLE]

The rate of change of the saturation $s$ is given by the conservation equation:

[TABLE]

where $g(s)$ is a nonlinear function. The constraint (2) still holds for an incompressible fluid. Equation (32) is similar to Darcy’s law for single-phase flow. An adaptation of the methodology outlined in this paper should be suited to the solution of (32). The hyperbolic saturation equation (34) poses more problems, in particular because its nonlinear nature leads to the formation of a shock in the solution buckley . Previous studies king indicate how a path integral formulation of this equation can be formed, providing an opportunity for future research.

Acknowledgements.

MJEW was supported through a Janet Watson scholarship from the Department of Earth Science and Engineering and a studentship in the Centre for Doctoral Training on Theory and Simulation of Materials funded by the EPSRC (EP/L015579/1), both at Imperial College London. We acknowledge the support of the Imperial College Research Computing Service HPC . We thank Dr. Stephan Dürr for a helpful discussion.