Flow in bounded and unbounded pore networks with different connectivity
Daniel W. Meyer Artur Gomolinski
physics.flu-dyn nlin.AO 76S05 (Primary) 45H05 (Secondary)




#1 1






Flow in Bounded and Unbounded Pore Networks with Different Connectivity
    Daniel W. MeyerThe first author gratefully acknowledges Branko Bijeljic from the Imperial College London for providing network data. Moreover, he acknowledges financial support from ETH Zürich.
 and Artur Gomolinski

  Institute of Fluid Dynamics, ETH Zürich
    June 28, 2019
======================================================================================================================================================================================================================================================================




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.



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









§ 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 <cit.>, we have developed the non-local Darcy formulation <cit.>. More recently, <cit.> 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 <ref>), the following expression for the averaged pressure p(x) was derived

    p(x)∫_D T(x,x^')dx^' - ∫_D T(x,x^')p(x^') dx^' = q(x)

<cit.>. In expression eqnlD, 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.[Slabs are considered as regions of statistical space-stationarity with respect to flow and pressure.] More specifically based on figure <ref>, p(x_A) ≡ n_A^-1∑_I∈ 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) ≡ (C h)^-1∑_I∈ 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^') [L/(F T) = T/M] is defined as

    T(x_A,x_B) ≡1/C h^2∑_I∈ A∑_J∈ B T_I,J(p_I - p_J)/p(x_A) - p(x_B),

where I and J are indices of pores that reside in network slabs A and B, respectively, both having a thickness h (see figure <ref>). T_I,J [L^5/(F T)] is the hydraulic conductivity of the throat that connects pores I and J.[Since 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

    T_I,J≡π D_I,J^4/128 μ L_I,J,

where D_I,J and L_I,J are the diameter and length, respectively, of the throat connecting pores I and J. μ [M/(LT)] is the dynamic viscosity of the fluid in the network.) In the non-local Darcy formulation eqnlD, T(x,x^') represents the non-local pressure transmission effect stated earlier. In the definitions of p(x), q(x), and T(x,x^'), 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^') directly reflects the connectivity of a porous medium, which is an aspect that has been investigated in the past by different groups. For example, <cit.> 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 <cit.>. In the followup contribution <cit.>, their methodology is combined with the multirate mass transfer model (MRMT) <cit.> and model predictions are compared against reference simulations. Moreover, in the context of discrete fracture networks (DFN) <cit.>, <cit.> 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

    -k/μ^2 p(x)x^2 = q(x)  k ≡μ∫_0^∞ s^2 T(s)  s

[L^2] and the dynamic viscosity μ [M/(LT)] of the fluid in the network <cit.>.

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^', but solely on their separation distance s ≡ x - x^'. 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 eqT).

In our earlier work <cit.>, 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 <cit.>). 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 <ref> 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 eqT 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 <cit.>. The stationary setup will be discussed in detail and assessed with a validation study, where we compare the network permeability resulting via expression eqDarcyLimit 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 <ref>. The theory for bounded networks is presented and assessed in section <ref>.




§ SPACE-STATIONARY PORE NETWORKS AND FLOW




 §.§ 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 <ref>, 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 <cit.>. This algorithm accounts for correlation between throat thickness and size of connected pores as proposed earlier by <cit.>. 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 <cit.>. 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]× [0,L_y]× [0,L_z] for the periodic network to be generated is prescribed.[To 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.]

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

  * To establish periodic throat connections reaching over boundaries of the target volume, pores located in the slices 0 ≤ x < L_m and L_x-L_m ≤ x < L_x are copied/shifted in x-direction outside the target volume by +L_x and -L_x, respectively (see figure <ref>). 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.

  * 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,… 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.

  * Throat connections involving a pore copy (resulting from step <ref>) 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 <ref> are removed.

  * 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 <ref> 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.





 §.§ 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 <cit.>.) The corresponding development is based on the flux balance of an arbitrary pore I:

    ∑_J Q_I,J = 0  ∀ I  Q_I,J≡ T_I,J (p_I - p_J)

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^' and a mean part ⟨ p(x_I)⟩, i.e., p_I = p_I^' + ⟨ p(x_I)⟩, where the latter is linearly varying in x-direction, i.e.,

    ⟨ p(x)⟩ = P(1-x/L_x),

and determined by the pressure magnitude P. Insertion of this decomposition into flux balance eqQ leads to

    ∑_J T_I,J (p_I^' - p_J^') = P/L_x∑_J T_I,J (x_I - x_J)  ∀ I,

which is an equation system for the pressure fluctuations p_I^'. In order to fix the pressure level of the fluctuations, equation system eqps has to be modified for example by replacing the flux balance of pore I = 1 with the auxiliary condition p_1^' = 0. Note that the throat fluxes Q_I,J are determined based on total pore pressures and therefore are proportional to P.




 §.§ Extraction of Conductivity Distribution




  §.§.§ Formulation


From the network and flow topology resulting from the techniques outlined in the previous two sections, the conductivity distribution eqT can be extracted by decomposing the network into slabs of thickness h as sketched in figure <ref>. However, since both the network and the flow topology are spatially independent, for example T(s ≡ 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 ≡ x_A-x_B. A corresponding adaptation of expression eqT leads to

    T(s) = 1/C h^2h/L_x∑_i = 1^L_x/h ∑_i-1 ≤ x_I/h < i  ∑_i-1 ≤ (x_J-s)/h < i T_I,J(p_I - p_J)/s P/L_x,

with separations s being even multiples of h. In expression eqTs, the slab pressure difference from the denominator of equation eqT was determined based on the prescribed mean pressure drop P/L_x. Since the pore pressures p_I and p_J in equation eqTs are proportional to the magnitude P (see expressions eqpm and eqps), the effect of P on T(s) cancels.




  §.§.§ Convergence



Based on the sandstone pore-network data-set available in <cit.> and the algorithm from section <ref>, a periodic pore network of dimensions (14)^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 <ref> with a dynamic viscosity μ = 8.9× 10^-4 and conductivity distributions were extracted by using expression eqTs with different slab thicknesses h = L_m/n. The resulting distributions for n = 8, 32, and 128 are depicted in figure <ref>. 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× 10^6 h/L_x ≈ 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 eqDarcyLimit numerically approximated with a trapezoidal rule, the permeability k = 3.86× 10^-13^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

    k = μQ_x/L_y L_zL_x/P.

This exact agreement confirms the Darcy limit eqDarcyLimit implied by the non-local Darcy theory.




  §.§.§ 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 <cit.>, which induces a non-stationary flow field (see figure <ref>). This comparison will reveal two inconsistencies. Like in our earlier contribution <cit.>, a pore network slice of extensions 3× 30× 30 was generated using the previously outlined algorithm and the sandstone dataset from <cit.>. The resulting network is comprised of 1.34 million pores and 2.85 million throats.[The maximum throat length L_m remains unchanged. Furthermore, L_x/L_m ≈ 6.] Firstly, we apply a mean pressure gradient in line with section <ref> and determine T(s) based on expression eqTs. The resulting distribution is depicted in figure <ref> (blue dots) and is, despite the different network extensions, in agreement with the one from the cubical network discussed in the previous section <ref> (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 ≈ 1.5L_m (see left panel in figure <ref>). 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 <ref> 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× 10^-13^2, that results via expression eqDarcyLimit from T(s), deviates from k = 3.86× 10^-13^2 reported in connection with the space-stationary setup. Moreover, k = 4.64× 10^-13^2 is inconsistent with k = 4.19× 10^-13^2, where the latter was calculated with equation eqk 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 <ref> and <ref>, 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 <ref>. The newly advocated setup, on the other hand, displays a constant mean pressure drop in the whole network (figure <ref> 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.




§ 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 <ref>, 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):

    Q_A = ∑_I∈ A∑_J∈ R_1∪ S∪ R_2 T_I,J(p_I - p_J)
           = ∑_I∈ A[∑_B∈ S∑_J∈ B T_I,J(p_I - p_J) + ∑_D_1∈ R_1∑_J∈ D_1x_J-x_I/x_S1-x_I T_I,J(p_I - p_R1) . 
          .  + ∑_D_2∈ R_2∑_J∈ D_2x_J-x_I/x_S2-x_IT_I,J(p_I - p_R2)] 
           = ∑_B∈ S∑_I∈ A∑_J∈ B T_I,J(p_I - p_J) + ∑_D_1∈ R_1∑_I∈ A∑_J∈ D_1x_J-x_I/x_S1-x_I_≈x_D1-x_A/x_S1-x_A T_I,J(p_I_≈ p(x_A) - p_R1) + …
           = Ch ∑_B∈ S T(x_A,x_B)[p(x_A)-p(x_B)]h + ∑_D_1∈ R_1x_D1-x_A/x_S1-x_A∑_I∈ A∑_J∈ D_1 T_I,J_≈ Ch^2 T(x_A,x_D1)[p(x_A) - p_R1] + …
           = C[∫_x_S1^x_S2 T(x_A,x)[p(x_A)-p(x)]  x + [p(x_A) - p_R1] ∫_-∞^x_S1x-x_A/x_S1-x_AT(x_A,x)  x. 
          .  +  [p(x_A) - p_R2] ∫_x_S2^∞x-x_A/x_S2-x_AT(x_A,x)  x] h = 0.

In the first step of the derivation of flux Q_A, the domain R_1∪ S ∪ 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 eqTij for cylindrical throats. In the remaining steps, the definition of the conductivity distribution eqT 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≈ p(x_A). (The validity of this assumption is assessed in section <ref>.) In the last step of derivation eqQA, the Riemann sums over slabs in regions S, R_1, and R_2 were replaced for slab thicknesses h → 0 by integrals <cit.>. Summation of all fluxes Q_A over all slabs A in sample S leads to

    ∑_A∈ SQ_A/C =       
    ∫∫_x_S1^x_S2 T(x^',x)[p(x^')-p(x)]  x x^'   +   ∫_x_S1^x_S2 [p(x^') - p_R1] ∫_-∞^x_S1x-x^'/x_S1-x^'T(x^',x)  x  x^'
       +   ∫_x_S1^x_S2 [p(x^') - p_R2] ∫_x_S2^∞x-x^'/x_S2-x^'T(x^',x)  x  x^' = 0.

With T(x,x^') = T(x^',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^' 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.




 §.§ Mean Pore-Pressure Drop in the Sample


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

    p(x) = p_R1 + x-x_S1/L(p_R2-p_R1)

is a solution of equation eqQS. Insertion into equation eqQS leads to

    ∫_x_S1^x_S2 (x^' - x_S1) ∫_-∞^x_S1x-x^'/x_S1-x^'T(x^',x)  x  x^'   +   ∫_x_S1^x_S2 (x^' - x_S2) ∫_x_S2^∞x-x^'/x_S2-x^'T(x^',x)  x  x^'
    
     = ∫_x_S1^x_S2[∫_-∞^x_S1 (x-x^')T(x^',x)  x .    +   . ∫_x_S2^∞ (x-x^')T(x^',x)  x ]  x^'
    
     = ∫_x_S1^x_S2[∫_-∞^∞ (x-x^')T(x^',x)  x .    -   . ∫_x_S1^x_S2 (x-x^')T(x^',x)  x ]  x^',

where the first inner integral is equal to zero since it can be written as ∫_-∞^∞ x^'' T(x^',x^'+x^'')  x^'' = ∫_-∞^∞ x^'' T(|x^''|)  x^'' for a space-stationary network with T(x,x^') = T(|x-x^'|). Moreover, the second integral is equal to zero based on an analogous argument made after equation eqQS and thus the linear pressure drop eqpx indeed satisfies equation eqQS under the condition of a space-stationary conductivity distribution.




 §.§ 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 eqQS and written here with a space-stationary conductivity distribution,

    q_S,R1    =    ∫_x_S1^x_S2 [p(x^') - p_R1] ∫_-∞^x_S1x-x^'/x_S1-x^'T(x-x^')  x  x^'
        =    ∫_x_S1^x_S2p(x^') - p_R1/x_S1-x^'∫_-∞^x_S1-x^' s  T(s)  s  x^'
        =    ∫_0^L p(s^'+x_S1) - p_R1/-s^'∫_-∞^-s^' s  T(s)  s  s^'
        =    ∫_0^L [p(s^'+x_S1) - p_R1]1/s^'∫_s^'^∞ s  T(s)  s  s^'.

With the linear mean pressure drop eqpx we obtain

    q_S,R1 = p_R2-p_R1/L∫_0^L ∫_s^'^∞ s  T(s)  s  s^'.

After switching the integration order based on figure <ref>a, the double integral in expression eqqSR1 can be rewritten as

    ∫_0^L ∫_s^'^∞ s  T(s)  s  s^'   =   ∫_0^∞∫_0^min(s,L) s  T(s)  s^' s = ∫_0^∞min(s,L) s  T(s)  s 
       =   ∫_0^L s^2 T(s)  s + ∫_L^∞ Ls T(s)  s,

which will be useful later.




 §.§ Short-Circuit Flux between Reservoirs


So far, the derivations in section <ref> 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 <ref>, 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 eqQA, the short-circuit flux is estimated in the following. The flux between two slabs D_1 and D_2 is given as 

    Q_D1,D2 = C x_D2-x_D1/L T(x_D1-x_D2) (p_R1-p_R2) h^2

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

    q_R1,R2    =    1/C∑_D_1∈ R_1∑_D_2∈ R_2 Q_D1,D2
        =    p_R1-p_R2/L∫_-∞^x_S1∫_x_S2^∞ (x-x^') T(x-x^')  x   x^'
        =    p_R1-p_R2/L∫_-∞^x_S1∫_x_S2-x^'^∞ s T(s)  s  x^'
        =    p_R1-p_R2/L∫_L^∞∫_x_S2-s^x_S1 s T(s)  x^' s 
        =    p_R1-p_R2/L∫_L^∞ (s-L) s T(s)  s.

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

Finally, summing the two flux components, with q_R1,S≡ -q_S,R1 determined by equations eqqSR1 and eqqSR12, leads to

    q_R1,S + q_R1,R2    =    p_R1-p_R2/L[∫_0^L s^2 T(s)  s + ∫_L^∞ Ls T(s)  s + ∫_L^∞ (s-L) s T(s)  s ] 
        =    p_R1-p_R2/L∫_0^∞ s^2 T(s)  s = p_R1-p_R2/Lk/μ,

where the last equality results from the Darcy limit eqDarcyLimit. 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.




 §.§ 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 <ref>, we start our verification study with the same network type. By using the algorithm outlined in section <ref> and the sandstone dataset from <cit.>, a periodic network of extensions 5L_m× 59L_m× 59L_m with L_m = 0.489, 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 <ref>. The resulting pore pressures along the network are reported in figure <ref>. While the mean pressure (white line) follows a linear trend in good agreement with the theoretical result reported in section <ref>, local pore pressures vary significantly and are mostly far from the mean pressure.




  §.§.§ Exact Short-Circuit Flux




In order to assess the flux predictions from sections <ref> and <ref>, 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 <ref>. While the theoretical q_R1,S correctly transitions from 0 to the asymptotic value k/μ = 3.86× 10^-13^2/[8.9× 10^-4] = 4.34× 10^-9^3 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→ 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 <ref>). Consequently, the flux through the sample is simply given by

    q_R1,R2(L→ 0) = 1/C∑_x_I < 0∑_x_J > Lx_J - x_I/LT_I,J  (p_R1-p_R2) = p_R1-p_R2/C L∑_x_I < 0∑_x_J > L (x_J - x_I)T_I,J.

Therefore as L→ 0, the limiting flux eqqR1R20 becomes independent of the network thickness L for a constant pressure drop (p_R1-p_R2)/L.[In expression eqqR1R20, x_I and x_J correspond to the pore positions of the uncut throats that go through sample S.] In figure <ref>a, the limiting flux is included with a black circle being in agreement with the numerical limit seen in figure <ref>b. Motivated by the reasoning behind expression eqqR1R20, the approximate expression eqqR1R2 for q_R1,R2 can be replaced by an exact relation: We write instead of equation eqQD1D2

    Q^'_D1,D2 = ∑_I∈ D_1∑_J∈ D_2x_J-x_I^≈ x_D2-x_D1/L T_I,J(p_I-p_J_=p_R1-p_R2) = C x_D2-x_D1/L T^'(x_D1-x_D2) (p_R1-p_R2) h^2

with the geometric conductivity distribution defined as

    T^'(x_A-x_B) ≡1/Ch^2∑_I∈ A∑_J∈ B T_I,J.

T^'(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 <ref>a. Insertion of expression eqQtD1D2 in derivation eqqR1R2 in place of Q_D1,D2 leads to

    q^'_R1,R2 = p_R1-p_R2/L∫_L^∞ (s-L) s T^'(s)  s.

As seen in figure <ref>, this flux is in agreement with its numerical counterpart extracted from the sandstone network. Moreover, it is consistent with the flux limit eqqR1R20.




  §.§.§ 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 eqQA). In figure <ref>, 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 <cit.> displaying a higher pore connectivity than the sandstone. From the network data obtained from <cit.>, a beakpack network of extensions 4L_m× 55L_m× 55L_m with L_m = 1.96× 10^-4, 3× 10^5 pores, and one million throats was generated by following the steps listed at the beginning of section <ref>. 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× 41L_m× 41L_m with L_m = 4.89× 10^-4 and is comprised of 2× 10^5 pores and one million throats.





The pore pressure variability for the beadpack and the homogeneous network are provided in figures <ref> and <ref>, 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 <ref>. Accordingly, the theoretical results are in better agreement with the corresponding pore network data as documented in figures <ref> and <ref>. Even though the theoretical results in figure <ref>a 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 ≈ 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^'(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 eqqSR1 for q_R1,S (with q_R1,S = -q_S,R1) is given by

    q^'_R1,S≡p_R1-p_R2/L[ ∫_0^L s^2 T(s)  s + ∫_L^∞ Ls T^'(s)  s ],

which is included in figures <ref>a, <ref>a, and <ref>a (light red solid lines). Here, in the second term T(s) was replaced by T^'(s). This term is large for small L and vanishes as L→ L_m, while the first term going to k/μ 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^'(s) in expression eqqtR1S. Such long-range pressure effects are particularly pronounced in less-connected networks (see sandstone and beadpack results in figures <ref> and <ref>).




§ 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 slab[Slabs 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.


plainnat


