Flow in bounded and unbounded pore networks with different connectivity

TL;DR

This paper develops a method to extract hydraulic conductivity distributions and presents a theory for mean pressure and flow in bounded pore networks, validated against numerical results, with improvements for networks of varying connectivity.

Contribution

It introduces a robust method for hydraulic conductivity extraction and a new theory for flow in bounded networks, extending previous non-local Darcy law formulations.

Findings

The method accurately extracts hydraulic conductivities.

The theory predicts flow well in highly connected networks.

Improvements enhance accuracy for less connected networks.

Abstract

This work is concerned with the intricate interplay between node or pore pressures and connection or throat conductivities in flow or pore networks. A setting similar to pore networks is given by fracture networks. Recently, a non-local generalization of Darcy's law for flow and transport in porous media was presented in the context of unbounded or periodic pore networks. In this work, we first outline a robust method for the extraction of the hydraulic conductivity distribution, which is at the heart of the non-local Darcy formulation. Second, a theory for mean pressure and flow in bounded networks is outlined. Predictions of that theory are validated against numerical network results and it is demonstrated that the theory works well for networks with high connectivity involving pores with high coordination numbers. For other networks, improvements to the outlined theory are proposed and their accuracy is assessed.

Equations

$$p(x){\int }_{D}T(x,x^{\prime })d{x}^{\prime }-{\int }_{D}T(x,x^{\prime })p(x^{\prime })d{x}^{\prime }=q(x)$$

$$T(x_A,x_B)\equiv \frac{1}{C{h}^2}\frac{{\sum }_{I\in A}{\sum }_{J\in B}{T}_{I,J}(p_I-p_J)}{p(x_A)-p(x_B)},$$

$$T_{I,J}\equiv \frac{\pi D_{I,J}^4}{128\mu L_{I,J}},$$

$$-\frac{k}{\mu }\frac{{\mathrm{d}}^{2}p(x)}{\mathrm{d}{x}^2}=q(x)\text{mbox}withthepermeabilityk\equiv \mu {\int }_0^{\infty }{s}^{2}T(s)\mathrm{d}s$$

$$\sum _J{Q}_{I,J}=0\forall I\text{mbox}with{Q}_{I,J}\equiv T_{I,J}(p_I-p_J)$$

$$⟨p(x)⟩=P\left(1-\frac{x}{L_x}\right),$$

$$\sum _J{T}_{I,J}(p_I^{\prime }-p_J^{\prime })=\frac{P}{L_x}\sum _J{T}_{I,J}(x_I-x_J)\forall I,$$

$$T(s)=\frac{1}{C{h}^2}\frac{h}{L_x}\frac{\sum _{i=1}^{L_x/h}\sum _{i-1\le x_I/h<i}\sum _{i-1\le (x_J-s)/h<i}{T}_{I,J}(p_I-p_J)}{sP/L_x},$$

$$k=\mu \frac{Q_x}{L_y{L}_z}\frac{L_x}{P}.$$

$$Q_A=\sum _{I\in A}\sum _{J\in R_1\cup S\cup R_2}{T}_{I,J}(p_I-p_J)$$

$$\sum _{A\in S}\frac{Q_A}{C}=$$

$$\int {\int }_{x_{S1}}^{x_{S2}}T(x^{\prime },x)[p(x^{\prime })-p(x)]\mathrm{d}x\mathrm{d}{x}^{\prime }$$

$$p(x)=p_{R1}+\frac{x-x_{S1}}{L}(p_{R2}-p_{R1})$$

$${\int }_{x_{S1}}^{x_{S2}}(x^{\prime }-x_{S1}){\int }_{-\infty }^{x_{S1}}\frac{x-x^{\prime }}{x_{S1}-x^{\prime }}T(x^{\prime },x)\mathrm{d}x\mathrm{d}{x}^{\prime }$$

$$={\int }_{x_{S1}}^{x_{S2}}[{\int }_{-\infty }^{x_{S1}}(x-x^{\prime })T(x^{\prime },x)\mathrm{d}x$$

$$={\int }_{x_{S1}}^{x_{S2}}[{\int }_{-\infty }^{\infty }(x-x^{\prime })T(x^{\prime },x)\mathrm{d}x$$

$$q_{S,R1}$$

$$q_{S,R1}=\frac{p_{R2}-p_{R1}}{L}{\int }_0^L{\int }_{s^{\prime }}^{\infty }sT(s)\mathrm{d}s\mathrm{d}{s}^{\prime }.$$

$${\int }_0^L{\int }_{s^{\prime }}^{\infty }sT(s)\mathrm{d}s\mathrm{d}{s}^{\prime }$$

$$Q_{D1,D2}=C\frac{x_{D2}-x_{D1}}{L}T(x_{D1}-x_{D2})(p_{R1}-p_{R2})h^2$$

$$q_{R1,R2}$$

$$q_{R1,S}+q_{R1,R2}$$

$$q_{R1,R2}(L\to 0)=\frac{1}{C}\sum _{x_I<0}\sum _{x_J>L}\frac{x_J-x_I}{L}{T}_{I,J}(p_{R1}-p_{R2})=\frac{p_{R1}-p_{R2}}{CL}\sum _{x_I<0}\sum _{x_J>L}(x_J-x_I)T_{I,J}.$$

$$Q_{D1,D2}^{\prime }=\sum _{I\in D_1}\sum _{J\in D_2}\frac{\stackrel{\approx x_{D2}-x_{D1}}{\overbrace{x_J-x_I}}}{L}{T}_{I,J}(\underset{=p_{R1}-p_{R2}}{\underbrace{p_I-p_J}})=C\frac{x_{D2}-x_{D1}}{L}{T}^{\prime }(x_{D1}-x_{D2})(p_{R1}-p_{R2})h^2$$

$$T^{\prime }(x_A-x_B)\equiv \frac{1}{C{h}^2}\sum _{I\in A}\sum _{J\in B}{T}_{I,J}.$$

$$q_{R1,R2}^{\prime }=\frac{p_{R1}-p_{R2}}{L}{\int }_L^{\infty }(s-L)s{T}^{\prime }(s)\mathrm{d}s.$$

$$q_{R1,S}^{\prime }\equiv \frac{p_{R1}-p_{R2}}{L}\left[{\int }_0^L{s}^{2}T(s)\mathrm{d}s+{\int }_L^{\infty }Ls{T}^{\prime }(s)\mathrm{d}s\right],$$

Full text

Keywords: non-local Darcy, pore network, connectivity, conductivity distribution, boundary

1 Introduction

Based on a porous-medium representation that idealizes the pore space in the form of a set of pores connected through a set of throats [e.g., 4, 2], we have developed the non-local Darcy formulation [5]. More recently, Jenny and Meyer [9] have related this non-local formulation to a particle-based transport description, which generalizes existing continuous time random walk (CTRW) models. Unlike Darcy’s law, where the flow at one point is determined by the local pressure gradient, the non-local formulation accounts for short- as well as long-range pressure transmission due to throats connecting pores over a range of distances. More specifically, for three-dimensional flow-networks that are statistically homogeneous in the spatial directions normal to the mean flow in $x$-direction (see figure 1), the following expression for the averaged pressure $p(x)$ was derived

[TABLE]

[see equation (5) in 5]. In expression (1), $p(x)$ is the average pore pressure with unit $[F/L^{2}]$ resulting from pores in a thin slab of thickness $h$ centered at $x$.111Slabs are considered as regions of statistical space-stationarity with respect to flow and pressure. More specifically based on figure 1, $p(x_{A})\equiv n_{A}^{-1}\sum_{I\in A}p_{I}$, where $n_{A}$ is the number of pores in slab $A$, $p_{I}$ is the pressure in pore $I$, and index $I$ enumerates the subset of pores contained in slab $A$. Similarly the specific source term $q(x)$ $[(L^{3}/T)/L^{3}=1/T]$ is defined as $q(x_{A})\equiv(Ch)^{-1}\sum_{I\in A}Q_{I}$ with $Q_{I}$ $[L^{3}/T]$ being the net-flow out of pore $I$ ($Q_{I}>0$) and $C$ being the cross-sectional area of the pore network perpendicular to $x$. The conductivity distribution $T(x,x^{\prime})$ $[L/(FT)=T/M]$ is defined as

[TABLE]

where $I$ and $J$ are indices of pores that reside in network slabs $A$ and $B$, respectively, both having a thickness $h$ (see figure 1). $T_{I,J}$ $[L^{5}/(FT)]$ is the hydraulic conductivity of the throat that connects pores $I$ and $J$.222Since $T_{I,J}=T_{J,I}$, it follows that $T(x_{A},x_{B})=T(x_{B},x_{A})$. (In the case of cylindrical throats, the throat conductivity $T_{I,J}$ is given as

[TABLE]

where $D_{I,J}$ and $L_{I,J}$ are the diameter and length, respectively, of the throat connecting pores $I$ and $J$. $\mu$ $[M/(LT)]$ is the dynamic viscosity of the fluid in the network.) In the non-local Darcy formulation (1), $T(x,x^{\prime})$ represents the non-local pressure transmission effect stated earlier. In the definitions of $p(x)$, $q(x)$, and $T(x,x^{\prime})$, the slab thickness $h$ is assumed to be sufficiently small, such that variations in these quantities are resolved and at the same time large enough such that statistical errors inherent to finite pore counts remain small.

The conductivity distribution $T(x,x^{\prime})$ directly reflects the connectivity of a porous medium, which is an aspect that has been investigated in the past by different groups. For example, Tyukhova et al. [16] outline a methodology to identify a network of channels of least resistance or preferential flow in a Darcy-flow context. These channels result from an iterative deformation process that is driven by local differences in the heterogeneous hydraulic conductivity [16, section 2.1]. In the followup contribution [15], their methodology is combined with the multirate mass transfer model (MRMT) [6] and model predictions are compared against reference simulations. Moreover, in the context of discrete fracture networks (DFN) [10, 13], Maillot et al. [11] compare commonly-used Poisson DFNs (comprised of statistically independent fractures) to kinematic fracture models (KFM), where the fracture pattern results from a growth process subject to fracture interaction. They find that the Poisson and KFM DFNs display quite different fracture connectivities leading to significant differences in flow and permeability.

Coming back to the non-local Darcy formulation, we could show in the limit of scales of interest being much larger than the longest throats in the network, that the formulation reduces to classical Darcy flow as

[TABLE]

$[L^{2}]$ and the dynamic viscosity $\mu$ $[M/(LT)]$ of the fluid in the network [5, Section 2.3].

For practical reasons, our past and current research efforts have been focusing on space-stationary settings, that is setups where $T$ does not depend on the actual points $x$ and $x^{\prime}$, but solely on their separation distance $s\equiv x-x^{\prime}$. This, however, does not only require the network topology to be space stationary, but the pressure statistics have to be spatially independent as well (see expression (2)).

In our earlier work [5], the conductivity distribution was extracted based on a network, which was periodic in mean-flow $x$-direction. A mean flow was induced in the network by introducing designated in-/outflow slabs, respectively (see [5, figure 4]). In the indicated slabs the pore pressures were set to a constant high/low value respectively, which induced a driving pressure gradient in the remaining regions of the network (see figure 4 in this work). The pressure field, however, is different in the in-/outflow slabs (where equal pore pressures are prescribed) compared to the other regions of the network (variable pressures). For this reason, $T(s)$ was estimated via expression (2) from slabs that were separated from the in-/outflow slabs by buffer zones. The primary function of these buffer zones was the relaxation of inhomogeneities in flow and pressure. In this work, we present a space-stationary network setup similar to the one briefly mentioned in our earlier work [9]. The stationary setup will be discussed in detail and assessed with a validation study, where we compare the network permeability resulting via expression (4) from $T(s)$ against its counterpart calculable from global network quantities such as size, net flow, and pressure drop.

In addition to the space-stationary setup intended for the extraction of $T(s)$, a theory for non-local flow in samples with finite extensions is developed. This development implies that the flux through a finite pore network is composed of two components: the first component represents short-circuit throats that directly connect the in- and outflow boundaries without any intermediate pores. The second flux component is comprised of throats that connect the two boundaries through intermediate pores in the network. Both fluxes are calculable as a function of the network thickness $L$ via integral expressions from the conductivity distribution $T(s)$. We validate predictions of our theory against numerical results stemming from networks of different topology.

After having outlined the non-local Darcy formulation in this introductory section, we present and validate the methodology for the construction of space-stationary pore networks and flow solutions in section 2. The theory for bounded networks is presented and assessed in section 3.

2 Space-Stationary Pore Networks and Flow

2.1 Space-Stationary Pore Networks

In order to imitate an unbounded space-stationary pore network, we resort to networks that are periodic in all spatial directions. For example in figure 1, a network with periodic throats in $x$-direction is depicted. To generate such networks, we use a slightly generalized version of the algorithm outlined by Idowu [8, Figure 3.3]. This algorithm accounts for correlation between throat thickness and size of connected pores as proposed earlier by [e.g., 4]. More specifically, a periodic network is resulting from the following steps:

Given is a base network in a cube volume that may stem from a sample of a natural porous medium [e.g., 14, Appendix I]. The base network is composed of pores that are connected through throats. Each pore has a position, radius, and list of connecting throats. Each throat is determined based on a radius and a pair of connecting pores. Moreover, a target volume $[0,L_{x}]\times[0,L_{y}]\times[0,L_{z}]$ for the periodic network to be generated is prescribed.333To avoid in steps 4 to 6 the formation of throats that connect a pore with itself, the side lengths $L_{x}$, $L_{y}$, and $L_{z}$ of the target volume need to be larger than the length $L_{m}$ of the longest throat in the base network. 2. 2.

With the same pore density as in the base network, pores are randomly placed in the target volume following a spatially uniform distribution. The properties of the newly placed pores such as radius and coordination number, i.e., the number of throats connected with a certain pore, are assigned by randomly drawing (with replacement) pores from the base network. 3. 3.

To establish periodic throat connections reaching over boundaries of the target volume, pores located in the slices $0\leq x<L_{m}$ and $L_{x}-L_{m}\leq x<L_{x}$ are copied/shifted in $x$-direction outside the target volume by $+L_{x}$ and $-L_{x}$, respectively (see figure 1). Here, $L_{m}$ is the length of the longest throat in the base network. The copy/shift operation is repeated in $y$- and subsequently in $z$-direction while including previously copied pores. 4. 4.

The first pore in the target volume with index $I=1$ and coordination number $c_{1}$ is connected through throats to the $c_{1}$ closest neighboring pores. This step is repeated for the other pores with indices $I=2,3,\ldots$ in the target volume while accounting for existing connections of neighboring pores. More specifically, a neighboring pore $J$, which has already been connected to $c_{J}$ pores or pore copies, is not connected to pore $I$. Instead the search neighborhood for connection candidates to pore $I$ is expanded up to a maximal distance $L_{m}$ from pore $I$. 5. 5.

Throat connections involving a pore copy (resulting from step 3) are rerouted to the original pore while storing the throat length that resulted from the location of the pore copy. With all throats connecting now pores in the target volume, the pore copies from step 3 are removed. 6. 6.

Throat radii are assigned such that throats connecting pores with large radii have a large radius as well. More specifically, for each throat connection that resulted from step 4 the radii of the connected pores are summed and a throat is randomly drawn from the base network. The radius of the throat with the largest pore-radius sum is set to the largest throat radius from all randomly drawn throats. Similarly, the radius of the throat with the second-largest pore-radius sum is set to the second-largest throat radius from the base network and so on for the remaining throat radii.

2.2 Space-Stationary Flow

Given the periodic network topology resulting from the algorithm outlined in the previous section, a space-stationary flow is induced based on a decomposition of the pressure into a sloping mean field and a periodic pressure fluctuation. (The same idea has been successfully applied in the context of Darcy flow in formations with heterogeneous permeability distributions [e.g., 3, 12].) The corresponding development is based on the flux balance of an arbitrary pore $I$:

[TABLE]

and where index $J$ enumerates the pores that are connected through a throat with conductivity $T_{I,J}$ to pore $I$. The pore pressure $p_{I}$ is now decomposed into a fluctuating part $p_{I}^{\prime}$ and a mean part $\langle p(x_{I})\rangle$, i.e., $p_{I}=p_{I}^{\prime}+\langle p(x_{I})\rangle$, where the latter is linearly varying in $x$-direction, i.e.,

[TABLE]

and determined by the pressure magnitude $P$. Insertion of this decomposition into flux balance (5) leads to

[TABLE]

which is an equation system for the pressure fluctuations $p_{I}^{\prime}$. In order to fix the pressure level of the fluctuations, equation system (7) has to be modified for example by replacing the flux balance of pore $I=1$ with the auxiliary condition $p_{1}^{\prime}=0$. Note that the throat fluxes $Q_{I,J}$ are determined based on total pore pressures and therefore are proportional to $P$.

2.3 Extraction of Conductivity Distribution

2.3.1 Formulation

From the network and flow topology resulting from the techniques outlined in the previous two sections, the conductivity distribution (2) can be extracted by decomposing the network into slabs of thickness $h$ as sketched in figure 1. However, since both the network and the flow topology are spatially independent, for example $T(s\equiv x_{A}-x_{B})=T(x_{A},x_{B})$ can be estimated not only from throats connecting two particular slabs centered at $x_{A}$ and $x_{B}$, but from all throats connecting $L_{x}/h$ pairs of slabs separated by $s\equiv x_{A}-x_{B}$. A corresponding adaptation of expression (2) leads to

[TABLE]

with separations $s$ being even multiples of $h$. In expression (8), the slab pressure difference from the denominator of equation (2) was determined based on the prescribed mean pressure drop $P/L_{x}$. Since the pore pressures $p_{I}$ and $p_{J}$ in equation (8) are proportional to the magnitude $P$ (see expressions (6) and (7)), the effect of $P$ on $T(s)$ cancels.

2.3.2 Convergence

Based on the sandstone pore-network data-set available in [7] and the algorithm from section 2.1, a periodic pore network of dimensions $(14\mbox{mm})^{3}$ with roughly 1.36 million pores and 2.9 million throats was generated. Next, a constant mean flow was induced using the setup outlined in section 2.2 with a dynamic viscosity $\mu=8.9\times 10^{-4}\mbox{kg/(ms)}$ and conductivity distributions were extracted by using expression (8) with different slab thicknesses $h=L_{m}/n$. The resulting distributions for $n=8$, 32, and 128 are depicted in figure 2. While $n=8$ slabs resolve the variability in $T(s)$ only poorly, 128 slabs lead to noisy estimates for $T(s)$ given on average $1.36\times 10^{6}h/L_{x}\approx 373$ pores per slab. With 32 slabs per maximal throat length $L_{m}$, the conductivity distribution is well resolved and statistical noise is kept at a low level. Based on that optimal $T(s)$ and expression (4) numerically approximated with a trapezoidal rule, the permeability $k=3.86\times 10^{-13}\mbox{m}^{2}$ is calculable. The same permeability is obtained from Darcy’s law with the total volumetric flux in $x$-direction $Q_{x}$ and the pressure gradient $P/L_{x}$ via

[TABLE]

This exact agreement confirms the Darcy limit (4) implied by the non-local Darcy theory.

2.3.3 Comparison with In-/Outflow Slab Method

In a next step, we compare the previously computed conductivity distribution against the one resulting from our earlier setup with in-/outflow slabs [5], which induces a non-stationary flow field (see figure 4). This comparison will reveal two inconsistencies. Like in our earlier contribution [5], a pore network slice of extensions $3\mbox{mm}\times 30\mbox{mm}\times 30\mbox{mm}$ was generated using the previously outlined algorithm and the sandstone dataset from Idowu [7]. The resulting network is comprised of 1.34 million pores and 2.85 million throats.444The maximum throat length $L_{m}$ remains unchanged. Furthermore, $L_{x}/L_{m}\approx 6$. Firstly, we apply a mean pressure gradient in line with section 2.2 and determine $T(s)$ based on expression (8). The resulting distribution is depicted in figure 3 (blue dots) and is, despite the different network extensions, in agreement with the one from the cubical network discussed in the previous section 2.3.2 (solid black line).

For the estimation of the conductivity distribution with the in-/outflow slab method, the in- and outflow slabs were embedded in buffer layers of width $L_{x}/4\approx 1.5L_{m}$ (see left panel in figure 4). Therefore, when estimating $T(s)$ only throats located outside of these layers in the remaining half of the network volume are accounted for. The resulting $T(s)$ is included in figure 3 through the red dashed line. It is seen that despite the buffer layers, inhomogeneities originating from the in-/outflow slabs induce a bias in the conductivity distribution $T(s)$ when comparing against the space-stationary setup (blue dots). As a result, the permeability $k=4.64\times 10^{-13}\mbox{m}^{2}$, that results via expression (4) from $T(s)$, deviates from $k=3.86\times 10^{-13}\mbox{m}^{2}$ reported in connection with the space-stationary setup. Moreover, $k=4.64\times 10^{-13}\mbox{m}^{2}$ is inconsistent with $k=4.19\times 10^{-13}\mbox{m}^{2}$, where the latter was calculated with equation (9) based on the net flux $Q_{x}$ in $x$-direction and with $P$ representing the pressure magnitude applied at the in-/outflow pores with positive/negative sign, respectively.

To further illustrate the difference between the space-stationary and the in-/outflow slab setups, pressure statistics are reported in figures 4 and 5, respectively. In the in-/outflow slabs located at $x=L_{x}/4$ and $3L_{x}/4$, respectively, constant pore pressures with zero variance were prescribed, which induce approximately linear pressure drops in large parts of the network, but at the same time inhomogeneities in the pore pressure statistics as shown in figure 4. The newly advocated setup, on the other hand, displays a constant mean pressure drop in the whole network (figure 5 left panel) and periodic pore-pressure fluctuations with no systematic dependence on the downstream coordinate $x$ (right panel). For networks of increasing size $L_{x}$ and buffer layer thickness $L_{x}/4$, it is expected that the flow field inhomogeneities in regions outside of the buffer layers will vanish and $T(s)$ from the in-/outflow method will converge to the distribution from the space-stationary setup.

3 Bounded Pore Networks

After having outlined a method for the unbiased estimation of the conductivity distribution in space-stationary settings, we turn in this section to a theoretical extension of the formulation for bounded networks. Based on the terminology from figure 6, we decompose the net flux $[L^{3}/T]$ from slab $A$ into a component going to other slabs in network sample $S$, and two components going to reservoirs $R_{1}$ and $R_{2}$ (pores cut away):

[TABLE]

In the first step of the derivation of flux $Q_{A}$, the domain $R_{1}\cup S\cup R_{2}$ was decomposed into slabs residing in sample $S$ and reservoirs $R_{1}$ and $R_{2}$. At the same time the conductivities $T_{I,J}$ of cut throats were corrected by their reductions in length consistent with expression (3) for cylindrical throats. In the remaining steps, the definition of the conductivity distribution (2) was used and it was assumed that the pore pressures $p_{I}$ and $p_{J}$ can be approximated with their respective slab mean pressures, e.g., $p_{I}\approx p(x_{A})$. (The validity of this assumption is assessed in section 3.4.) In the last step of derivation (3), the Riemann sums over slabs in regions $S$, $R_{1}$, and $R_{2}$ were replaced for slab thicknesses $h\to 0$ by integrals [5, equation (4)]. Summation of all fluxes $Q_{A}$ over all slabs $A$ in sample $S$ leads to

[TABLE]

With $T(x,x^{\prime})=T(x^{\prime},x)$, the first two-dimensional integral on the right hand side vanishes, since contributions from point pairs $(x_{1},x_{2})$ and $(x_{2},x_{1})$ in the quadratic $x$-$x^{\prime}$ integration region cancel as $T(x_{1},x_{2})[p(x_{1})-p(x_{2})]=-T(x_{1},x_{2})[p(x_{2})-p(x_{1})]$. The remaining two integrals represent the net fluxes from the sample $S$ into the reservoirs at pressures $p_{R1}$ and $p_{R2}$, respectively.

3.1 Mean Pore-Pressure Drop in the Sample

In a next step, we verify whether the linear mean pressure drop

[TABLE]

is a solution of equation (11). Insertion into equation (11) leads to

[TABLE]

where the first inner integral is equal to zero since it can be written as $\int_{-\infty}^{\infty}x^{\prime\prime}T(x^{\prime},x^{\prime}+x^{\prime\prime})\mathrm{d}x^{\prime\prime}=\int_{-\infty}^{\infty}x^{\prime\prime}T(|x^{\prime\prime}|)\mathrm{d}x^{\prime\prime}$ for a space-stationary network with $T(x,x^{\prime})=T(|x-x^{\prime}|)$. Moreover, the second integral is equal to zero based on an analogous argument made after equation (11) and thus the linear pressure drop (12) indeed satisfies equation (11) under the condition of a space-stationary conductivity distribution.

3.2 Flux from Sample to Reservoir

Based on the mean pressure drop found in the previous section, a compact expression for the flux per unit area $q_{S,R1}$ $[L^{3}/(TL^{2})]$ going from sample $S$ to, e.g., reservoir $R_{1}$ shall be derived: Starting point is the corresponding integral identified in equation (11) and written here with a space-stationary conductivity distribution,

[TABLE]

With the linear mean pressure drop (12) we obtain

[TABLE]

After switching the integration order based on figure 7a, the double integral in expression (14) can be rewritten as

[TABLE]

which will be useful later.

3.3 Short-Circuit Flux between Reservoirs

So far, the derivations in section 3 have focused on flow through sample $S$ via pores located in the sample. However, there is a second component to the flow between reservoirs $R_{1}$ and $R_{2}$ through sample $S$. As illustrated in figure 6, if $L<L_{m}$, there may be throats that go straight through sample $S$, short-circuiting reservoirs $R_{1}$ and $R_{2}$. Based on similar arguments as applied in the derivation of equation (3), the short-circuit flux is estimated in the following. The flux between two slabs $D_{1}$ and $D_{2}$ is given as

[TABLE]

and summation of all pairs of slabs $D_{1}$ and $D_{2}$ in reservoirs $R_{1}$ and $R_{2}$, respectively, leads to the net flux between $R_{1}$ and $R_{2}$ bypassing pores in samples $S$, i.e.,

[TABLE]

To get from the third to the fourth equality, the order of the integrations over the region depicted in figure 7b was switched.

Finally, summing the two flux components, with $q_{R1,S}\equiv-q_{S,R1}$ determined by equations (14) and (15), leads to

[TABLE]

where the last equality results from the Darcy limit (4). Interestingly, given the assumptions made, the total flux from reservoir $R_{1}$ through sample $S$ to reservoir $R_{2}$ is consistent with Darcy’s law.

3.4 Comparison with Pore-Network Results

In a last part, the theoretical results derived in the previous sections shall be verified against numerical results from pore network calculations. In line with the sandstone network considered in section 2.3, we start our verification study with the same network type. By using the algorithm outlined in section 2.1 and the sandstone dataset from Idowu [7], a periodic network of extensions $5L_{m}\times 59L_{m}\times 59L_{m}$ with $L_{m}=0.489\mbox{mm}$, one million pores, and 2.1 million throats was generated (corresponding to an average coordination number of 4.2). In a next step, periodic throats in $x$-direction were cut open at the network boundaries at $x=0$ and $L_{x}$ and constant pressures $p_{R1}$ and $p_{R2}$, respectively, were applied at the throat cutting points as sketched in figure 6. The resulting pore pressures along the network are reported in figure 8. While the mean pressure (white line) follows a linear trend in good agreement with the theoretical result reported in section 3.1, local pore pressures vary significantly and are mostly far from the mean pressure.

3.4.1 Exact Short-Circuit Flux

In order to assess the flux predictions from sections 3.2 and 3.3, the pore network is gradually cut back leading to networks of decreasing thickness $L_{x}$ or $L$. The resulting normalized fluxes are compared with their theoretical counterparts in figure 9. While the theoretical $q_{R1,S}$ correctly transitions from 0 to the asymptotic value $k/\mu=3.86\times 10^{-13}\mbox{m}^{2}/[8.9\times 10^{-4}\mbox{kg/(ms)}]=4.34\times 10^{-9}\mbox{m}^{3}\mbox{s/kg}$ consistent with Darcy’s law, it deviates from the numerical result at intermediate $L$. Moreover, $q_{R1,R2}$ is considerably higher in networks with $L\to 0$ compared to the theoretical prediction. For a network with vanishing thickness $L$, the likelihood to observe pores in the sample goes to zero (compare figure 6). Consequently, the flux through the sample is simply given by

[TABLE]

Therefore as $L\to 0$, the limiting flux (19) becomes independent of the network thickness $L$ for a constant pressure drop $(p_{R1}-p_{R2})/L$.555In expression (19), $x_{I}$ and $x_{J}$ correspond to the pore positions of the uncut throats that go through sample $S$. In figure 9a, the limiting flux is included with a black circle being in agreement with the numerical limit seen in figure 9b. Motivated by the reasoning behind expression (19), the approximate expression (17) for $q_{R1,R2}$ can be replaced by an exact relation: We write instead of equation (16)

[TABLE]

with the geometric conductivity distribution defined as

[TABLE]

$T^{\prime}(s)$ is determined by the network geometry, but in contrast to $T(s)$ it does not depend on the flow. A comparison of the flow- vs. the geometry-based conductivity distribution for the sandstone network is provided in figure 10a. Insertion of expression (20) in derivation (17) in place of $Q_{D1,D2}$ leads to

[TABLE]

As seen in figure 9, this flux is in agreement with its numerical counterpart extracted from the sandstone network. Moreover, it is consistent with the flux limit (19).

3.4.2 Pressure and Pore Connectivitiy

The remaining deviations in $q_{R1,S}$ are analyzed next by inspecting the simplifying assumption of setting pore pressures equal to mean slab pressures (compare derivation (3)). In figure 8, it was found for a sandstone network of thickness $L/L_{m}=5$ that pore pressures deviate significantly from the local mean pressure. A similar analysis as with the sandstone network was performed with a beadpack network [1] displaying a higher pore connectivity than the sandstone. From the network data obtained from Bijeljic [1], a beakpack network of extensions $4L_{m}\times 55L_{m}\times 55L_{m}$ with $L_{m}=1.96\times 10^{-4}\mbox{m}$, $3\times 10^{5}$ pores, and one million throats was generated by following the steps listed at the beginning of section 3.4. Moreover, a homogeneous network, where each pore is connected to ten neighboring pores through throats of equal radii, was studied. This network has extensions of $2L_{m}\times 41L_{m}\times 41L_{m}$ with $L_{m}=4.89\times 10^{-4}\mbox{m}$ and is comprised of $2\times 10^{5}$ pores and one million throats.

The pore pressure variability for the beadpack and the homogeneous network are provided in figures 11 and 13, respectively. These results illustrate that an increased coordination number (average values of 6.6 for the beadpack and 10 for the homogeneous network) leads to more homogeneous pore pressures in better agreement with the assumptions made in our theoretical analysis. Moreover, the geometry- and flow-based conductivity distributions are more alike than for the sandstone network as partly seen in figure 10. Accordingly, the theoretical results are in better agreement with the corresponding pore network data as documented in figures 12 and 14. Even though the theoretical results in figure 14a come close to their numerical counterparts in panel b of that figure, $q_{R1,S}$ from the network does not only approach the Darcy limit from below, as predicted by the theory, but for $L\approx L_{m}$ exceeds the Darcy limit. This is related to the fact that for thin networks, pore pressures are closely correlated with the reservoir pressures (conductivity is quantified by $T^{\prime}(s)$) and thus vary locally less than in the space-stationary section emerging for large $L$ in the center of the sample ($T(s)$ for conductivity). Thus, an ad-hoc improvement over the model (14) for $q_{R1,S}$ (with $q_{R1,S}=-q_{S,R1}$) is given by

[TABLE]

which is included in figures 9a, 12a, and 14a (light red solid lines). Here, in the second term $T(s)$ was replaced by $T^{\prime}(s)$. This term is large for small $L$ and vanishes as $L\to L_{m}$, while the first term going to $k/\mu$ becomes dominant. Despite this modification, pore pressures remain correlated over several throat generations in flow direction $x$ and thus the effect under inspection spans over lengths $L>L_{m}$, which is, however, not accounted for by $T^{\prime}(s)$ in expression (23). Such long-range pressure effects are particularly pronounced in less-connected networks (see sandstone and beadpack results in figures 9 and 12).

4 Summary

In the first part of this work, we have presented a method for the extraction of the flow-based conductivity distribution, which is at the heart of the non-local Darcy formulation. The conductivity distribution accounts for non-local pressure transmission resulting from throats which connect distant pores. Unlike the in-/outflow slab method, which induces a mean flow in the network through a singular source/sink pair, the new approach based on a space-stationary network and flow field does not require buffer layers that absorb or dampen flow non-stationarities resulting from the localized source/sink. A convergence study for different slab666Slabs are pore network sub-volumes with flow and pressure statistics that are to a good approximation statistically homogeneous in space. widths was successfully conducted and consistency with the classical Darcy limit has been verified.

By assuming a stationary conductivity distribution and by approximating pore pressures based on their local slab means, a theory describing the flow in bounded networks has been derived. The theory implies a constant mean pressure gradient in the sample and a net flux through the sample both consistent with Darcy’s law. In a next step, these theoretical results were verified with numerical network simulations. It was found that the stated assumptions are reasonable for highly-connected networks, but inadequate with realistic natural pore networks leading to inaccurate predictions. While approximately linear mean pore-pressure drops in agreement with the theory were found, high local pore-pressure variances and deviations in the total flux through thin networks were observed. In natural networks, the stated pressure variations are of similar size as the mean pressure drop over a characteristic network length (e.g., length of the longest throat $L_{m}$). Since the pressure is fixed to a constant value at boundaries (inducing zero variance), the flow statistics near the bounds deviate from the interior of the network, where the flow is similar to a periodic network (with non-zero variance). Introduction of a second conductivity distribution based entirely on the network geometry, and thus accounting for the zero pressure variability at boundaries, leads to reasonably accurate theoretical flux predictions. In conclusion, the presented theory is expected to provide accurate results for networks where the conductivity distribution based on flow is similar to the one based on the network geometry, with the latter ignoring the coupling between pore pressures and throat conductivities.

One remaining open aspect of this work is the organization of the flow along sequentially linked throats that connect the interior of the network with its boundaries. This effect is not captured by the present theory based on the two listed conductivity distributions, but clearly observable in the numerical network results: While in the theory, changes in the fluxes were predicted up to networks of maximal thickness given by the longest throat $L_{m}$ (both conductivity distributions have a spatial support of $2L_{m}$), in the numerical pore network results, flux dependencies for much thicker networks were found. It seems necessary to define an effective throat length that is longer than actual throats connecting just two pores. A simple example of the occurrence of an effective throat is implied by a throat pair that is linked through a common pore, which in turn has no further throat connections. Here, the effective throat results from the sequential combination of the throat pair, while the intermediate pore can be safely ignored.