Pressure and flow statistics of Darcy flow from simulated annealing

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.


I. INTRODUCTION
The study of flow in porous media has a wide range of applications, including hydrology (e.g. [1,2]), geothermal engineering (e.g. [3]), materials science (e.g. [1]), and the medical sciences (e.g. [4]). 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 [5]. 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") q of the fluid, and the pressure gradient ∇p: Here, K(x) = k(x)/µ, where k(x) is the effective permeability of the medium and µ 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 [6], Darcy's law can be recovered from the Navier-Stokes equations [7]. In this paper, we focus on single-phase, incompressible flow: In reservoir engineering, Darcy's law is solved numerically using the finite-volume method [8,9], while the finite-element method is common in hydrology.
In previous work [10], 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), from which we immediately obtain where q 0 is a constant. In higher dimensions, the equation corresponding to (3), has the general solution where f is a twice differentiable vector function and ∇ 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 [11,12]. 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 singlephase 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 [13][14][15][16][17]. 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 [18].
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.

A. Effective permeability
The effective permeability is modelled as a stochastic process. Various such models exist, with an overview is given in [19]. We have opted to model K(x) as a lognormal process, where L(x) is a Gaussian process with mean zero, characterized by its variance σ 2 and correlation length ξ. More information on this process can be found in [20]. This conventional choice [21] has the advantage of a strictly positive permeability. We emphasize, however, that our method applies to any choice of permeability distribution.

B. Path integral in higher dimensions
In an earlier paper [10], 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 deltafunctional: The term encodes the correlated Gaussian probability distribution of the log-permeability, where C L denotes the log-permeability covariance matrix.
The factor e − i Li is the Jacobian associated with integrating over the L i rather than over the e Li . Upon integration over the L i we obtain a path integral with the discrete "action" [10] Discrete pressure paths are generated according to the probability density (10). An analogous path integral to (8) in two dimensions is obtained with the standard procedure for classical statistical dynamics [22][23][24][25]: The delta-function enforces the discrete form of Darcy's law (5), where we have used the notation 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 [22][23][24][25] is to apply a functional Fourier transform, which yields a complex action. This is appropriate for formal studies in-volving 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: This expression is readily generalized to three dimensions.
Averages of pressure and correlation functions can be calculated from (14) by first generating permeability fields, then setting t to some value, and finally using the Metropolis-Hasting (MH) algorithm to minimize the discrete "action": 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.

C. Path integral for simulated annealing
To account for any number of dimensions, let us write the action in its continuum form: where the integral is carried out over the entire volume under consideration. Simulated annealing aims to find the pressure p(x) that minimizes the action (16) for fixed permeability K(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 δp and imposing the condition Any boundary conditions are unchanged, so δp(x) = 0 at the boundaries of the volume V . If there exists a p * such that the stationarity condition (17) holds, the action S is stationary at p * .
Retaining δp only to first order and performing an integration by parts, we obtain 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 which is Darcy's law (5) for incompressible flow. Thus, if p * can be found such that the stationarity condition is met at fixed K, then p * 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 [26].

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 [L 2 T −1 ], the flow q [L T −1 ] and the pressure (also known as the hydraulic head) p [L]. The total flow in a given direction is denoted by Our three-dimensional, rectangular prismatic porous medium is simulated on a grid of lattice elements, representing a domain of size The correlation lengths of the permeability field are 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 σ 2 , the defining feature of the six parameter sets, are given in Table I. The values in this table were those used by Nowak et al. [27], which enables a qualitative comparison to be made between our two approaches.

A. Simulation of the permeability fields
To simulate the permeability field, the first step is to generate the three-dimensional logpermeability field. Once the Gaussian field L(x) is available, definition (7) can be invoked to calculate the permeability field K(x). To generate L(x), we have made use of the circulant embedding technique [28], where the correlation matrix C of the desired field is embedded into a matrix  [30].

B. Simulated annealing
The simulated annealing algorithm is based on the MH algorithm, a step-by-step explanation of which can be found in [32]. 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.
1. Initialize a random pressure that is consistent with the boundary conditions.
2. Execute the MH algorithm some M 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(x)) and does not get stuck in a local minimum of the state space.
3. After every N s 1 steps, check the value of S. When the value of S starts to fluctuate around a constant value, go to step 4.

Adapt the MH accept/reject criterion to
"accept the change in the action with probability min(1, e −δS/T ) for some constant 0 < T < 1". This is a "cooling step" [31]. 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.
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 2 , say 2 = 10 −2 .
To expedite the simulated annealing algorithm, we made use of a technique known as overrelaxation (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 [33,34]. 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 lies "on the other side" of the minimum of the action: p * i is the value of p i that minimizes S, with all other parameters kept fixed. Since the overrelaxation 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 where k indicates the cooling step. The MH algorithm was executed M = 2, 000 times for all parameters. For σ 2 ≤ 1, repeating the cooling algorithm N s = 3, 000 times was found to be a good choice. For σ 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 floating point operations ("flops"). The key calculation in the FVM is a sparse matrix inversion. The sparse matrix solver UMFPACK [35] can solve such an equation in flops. Contrary to the FVM, the simulated annealing requires 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 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 [36][37][38] and directed sampling [39][40][41].

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 σ 2 . Based on information about boundedness, we have made parametric fits using the choices made in [27] for guidance.
All quantities were normalized for straightforward interpretation. The Dirichlet boundary conditions were chosen to yield a pressure difference The flow components were normalized as where K e is the theoretical expectation value of the permeability, I 0 = ∆p/L y and A = L x × 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 lognormal distribution is an obvious contender for parametric fits to the pressure. The probability density function is given by: 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 [42], as the generalized Central Limit Theorem predicts [43]. For a more extensive explanation of this theorem in the context of Darcy flow, see [20]. These boundary effects increase with the logpermeability variance, as can also be observed for the flow.   Table I. Squares represent the results of the finitevolume method, circles those of simulated annealing.
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 logpermeability variance, the flow statistics are most clearly log-normal. This is because the boundary effects are more strongly felt for high values of σ 2 , a pattern that can also be observed by comparing Figures 1 and 2 tion, 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 (σ /µ ) 2 → 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:  Figure 3. Normalized flow in main direction q * y for the log-permeability variances σ 2 given in Table I, measured at the point (0.5X, 0.5Y, 0.5Z). The log-normal pattern manifests itself most clearly for high values of σ 2 . For these high values, the boundary effects, which dictate the log-normality, have the strongest influence on the flow statistics.
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 [27]). 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.   Table I. The distribution of Q * y resembles that of q * y shown in Fig. 3. Due to the averaging over a cross-section, the distributions appear more Gaussian.
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 [45], or copulas, which describe the stochastic structure without reference to the corresponding marginal distributions [46].
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 remain-  Figure 5. Statistics of the normalized flow in the x-direction, q * x , again for the parameter set given in Table I. The parametric fits of the exponential power distribution reflect the symmetry of the statistics about zero. The tails become heavier for greater values of σ 2 . In the limit σ 1, the flow either continues along the main axis or "hits a wall" and reverses course. Thus, in this limit, the likelihood of small values for q * x is very small.
ing 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 finitevolume 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 twophase flow a generalized form of Darcy's law is used: 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: The rate of change of the saturation s is given by the conservation equation: 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 [47]. Previous studies [48] indicate how a path integral formulation of this equation can be formed, providing an opportunity for future research.