4.1. Model equations

To solve the model equations (2.9) and (2.10), the domain $\varOmega =[0,H_x]$ is discretized by an equidistant grid with $N$ cells of size $h_x=H_x/N$. A finite-volume method was employed where for each cell $\varOmega _I=[(I-1)h_x,Ih_x]$ ($I\in \{1,\ldots ,N\}$) the conditions

(4.1a)\begin{equation} A_{f_I}P_{f_I}=\sum_{J=1}^{N}\left\{G_{m_{IJ}}P_{m_J}\right\} +B_{l_I}P_l+B_{r_I}P_r \end{equation}

and

(4.1b)\begin{equation} A_{m_I}P_{m_I} =\sum_{J=1}^{N}\left\{G_{m_{JI}}P_{f_J}\right\} +C_{l_I}P_l +C_{r_I}P_r + K_m\left(P_{m_{I+1}}+P_{m_{I-1}}\right) \end{equation}

with

(4.2)\begin{gather} K_m=\frac{\bar{k}_m}{h_x^2}, \end{gather}

(4.3)\begin{gather}B_{l_I}=\frac{\hat{k}_f(x_I,x_l)}{x_I-x_l}, \end{gather}

(4.4)\begin{gather}B_{r_I}=\frac{\hat{k}_f(x_I,x_r)}{x_r-x_I}, \end{gather}

(4.5)\begin{gather}C_{l_I}=\int_{-\infty}^{x_l}{\hat{g}(x,x_I)}\,\textrm{d}x, \end{gather}

(4.6)\begin{gather}C_{r_I}=\int_{x_r}^{\infty}{\hat{g}(x,x_I)}\,\textrm{d}x, \end{gather}

(4.7)\begin{gather}G_{m_{IJ}}={\hat{g}(x_I,x_J)}h_x, \end{gather}

(4.8)\begin{gather}A_{f_I}=\sum_{J=1}^{N}\{G_{m_{IJ}}\} + B_{l_I} + B_{r_I}, \end{gather}

(4.9)\begin{gather}A_{m_I}=\sum_{J=1}^{N}\{G_{m_{JI}}\} + C_{l_I} + C_{r_I} + 2K_m, \end{gather}

(4.10a,b)\begin{gather}P_{m_{0}}=P_l\quad \textrm{and}\quad P_{m_{N+1}} = P_r \end{gather}

are enforced. The Gauss Seidel algorithm is used as a linear solver, i.e. one obtains the solution from the iteration equations

(4.11)\begin{equation} P^{\nu+1}_{f_I} =\frac{1}{A_{f_I}}\left( \sum_{J=1}^{I-1}\left\{ G_{m_{IJ}}P_{m_J}^{\nu+1} \right\} + \sum_{J=I}^{N}\left\{ G_{m_{IJ}}P_{m_J}^\nu \right\} + B_{l_I}P_l + B_{r_I}P_r\right) \end{equation}

and

(4.12)\begin{align} P^{\nu+1}_{m_I}&=\frac{1}{A_{m_I}}\left( \sum_{J=1}^{I}\left\{G_{m_{JI}}P_{f_J}^{\nu+1} \right\} + \sum_{J=I+1}^{N}\left\{ G_{m_{JI}}P_{f_J}^\nu \right\} + C_{l_I}P_l + C_{r_I}P_r \right.\nonumber\\ &\quad \left.\vphantom{\sum_{J=1}^{I}}+K_m\left(P_{m_{I-1}}^{\nu+1}+P_{m_{I+1}}^{\nu}\right)\right), \end{align}

where $P_l=p_b(0)$ and $P_r=p_b(H_x)$ are left and right boundary pressures. During each iteration $\nu +1$ the new pressure values $P_{m_I}^{\nu +1}$ and $P_{f_I}^{\nu +1}$ are computed in each cell in the order of the index $I$, and eventually, for $\nu \rightarrow \infty$, the solutions converge.

4.2. Darcy and continuity equations

To solve (2.5), the two-dimensional domain $\varOmega =[0,H_x]\times [0,H_y]$ is discretized by an equidistant, Cartesian grid with $N\times M$ cells $\varOmega _{I,J}=[(I-1)h_x,Ih_x]\times [(J-1)h_y,Jh_y]$, where $h_x=H_x/N$ and $h_y=H_y/M$. A finite-volume method is applied, which leads to the constraints

(4.13)\begin{equation} A_{I,J}P_{I,J}=A_{I-1,J}P_{I-1,J}+A_{I+1,J}P_{I+1,J}+A_{I,J-1}P_{I,J-1}+A_{I,J+1}P_{I,J+1}. \end{equation}

For the coefficients, harmonic averaging of the adjacent cell permeability values is employed, i.e. they read

(4.14)\begin{gather} A_{I-1,J}=\frac{2 K_{I,J}K_{I-1,J}}{K_{I,J}+K_{I-1,J}} \frac{1}{h_x^2}, \end{gather}

(4.15)\begin{gather}A_{I+1,J}=\frac{2 K_{I,J}K_{I+1,J}}{K_{I,J}+K_{I+1,J}} \frac{1}{h_x^2}, \end{gather}

(4.16)\begin{gather}A_{I,J-1}=\frac{2 K_{I,J}K_{I,J-1}}{K_{I,J}+K_{I,J-1}} \frac{1}{h_y^2}, \end{gather}

(4.17)\begin{gather}A_{I,J+1}=\frac{2 K_{I,J}K_{I,J+1}}{K_{I,J}+K_{I,J+1}} \frac{1}{h_y^2} \end{gather}

and

(4.18)\begin{equation} A_{I,J}=A_{I-1,J}+A_{I+1,J}+A_{I,J-1}+A_{I,J+1}. \end{equation}

Note that for all Monte Carlo studies presented in this paper fractures are represented as horizontal pixel lines with numerical permeabilities of $k^{numerical}_f=k_fa_f/h_y$, where $k_f$ and $a_f$ are fracture permeability and aperture, respectively. Periodic boundary conditions are applied in the direction of the second index (in the $y$-direction), i.e.

(4.19)\begin{equation} (P_{I,0},P_{I,M+1})=(P_{I,M},P_{I,1}) \end{equation}

and

(4.20)\begin{equation} (K_{I,0},K_{I,M+1})=(K_{I,M},K_{I,1}). \end{equation}

In the direction of the first index ($x$-direction), either Dirichlet boundary conditions are applied, i.e.

(4.21)\begin{equation} (P_{0,J},P_{N+1,J})=(2P_l-P_{1,J},2P_r-P_{N,J}) \end{equation}

and

(4.22)\begin{equation} (K_{0,J},K_{N+1,J})=(K_{1,J},K_{N,J}), \end{equation}

or periodic boundary conditions with an imposed mean pressure gradient of $-H_x^{-1}$, i.e.

(4.23)\begin{equation} (P_{0,J},P_{N+1,J})=(P_{N,J} + 1,P_{1,J} - 1) \end{equation}

and

(4.24)\begin{equation} (K_{0,J},K_{N+1,J})=(K_{N,J},K_{1,J}). \end{equation}

As for the model equations, the Gauss Seidel algorithm is used as a linear solver, i.e. the iteration equation

(4.25)\begin{equation} P^{\nu+1}_{I,J} =\frac{1}{A_{I,J}}\left( A_{I-1,J}P^{\nu+1}_{I-1,J}+A_{I+1,J}P^\nu_{I+1,J}+A_{I,J-1}P^{\nu+1}_{I,J-1}+A_{I,J+1}P^\nu_{I,J+1} \right) \end{equation}

is applied for all cells (of the order $I=1\rightarrow N$ and $J=1\rightarrow M$), until all pressure values $P_{I,J}^{\nu +1}$ converge.

5.1. Test case

Statistically homogeneous fractured rock samples are considered, whereas all quantities are assumed to be constant in the $z$-direction. Thus, for all the following Monte Carlo studies, the two-dimensional (2-D) domain $\varOmega =[0,H_x]\times [0,H_y]$ is used in place of the 3-D domain $\varOmega ^\textrm {3D}=[0,H_x]\times [0,H_y]\times [0,H_z=1]$. Further, all fractures are isolated, aligned with the mean flow (in the $x$-direction) and have the same length $l_f=0.25$. Stochasticity enters through the random fracture centre locations, i.e. for each Monte Carlo realization they were sampled from a uniform distribution with the specified density within a domain much larger than the one used for the computations. Fracture and matrix permeability are $k_f=2/a_f$ and $k_m=1$, respectively, where $a_f$ is the fracture aperture. The latter is assumed to be smaller than the grid resolution. Note that this implies that a single fracture is equally conductive as a pure matrix domain with $H_y=2$. Further, the mean fracture number density $\rho _f$ in the $x$–$y$-plane is $50$. An illustration of the test case is shown in figure 2. The fracture centres are marked with fuzzy dots, and from the right figure, it becomes apparent that also some fractures with their centres outside of the domain have to be taken into account. To generate independent samples, $10(H_x+l_f)H_y\rho _f$ fracture centres $\boldsymbol {x}_f$ are randomly placed (uniform distribution) in the sampling domain $\varOmega ^{sample}=[-H_x-l_f/2,H_x+l_f/2]\times [0,10H_y]$. Note that $\varOmega ^{sample}$ is much larger than $\varOmega$, which is to properly account for fluctuations of the local fracture density. If periodic boundary conditions are applied in the $x$-direction, all fractures with $\boldsymbol {x}_f\notin \varOmega$ are rejected, and for each accepted fracture a pixel line of length $l_f/h_x$ with permeability $k^{numerical}_f=k_fa_f/h_y$ is created, whereas periodicity is assumed for all pixels outside of $\varOmega$. If Dirichlet boundary conditions are applied, the accepted fractures are those with $\boldsymbol {x}_f\in [-l_f/2,H_x+l_f/2]\times [0,H_y]$; see figure 4. Regarding resolution and domain size, it was found that $h_x=h_y=0.01$ and $H_y=0.5$ are sufficient. Further, converged average flow rates with periodic boundary conditions (in the $x$-direction with a unit pressure drop) were achieved with $H_x=1.5$; convergence studies were performed for lengths up to $H_x=2.5$ and widths up to $H_y=1$. The applied grid resolution also ensures that the fractures of length $l_f=0.25$ can exactly be represented by lines of $25$ cells. Pressure surface plots of two realizations with $H_x=1.5$ and $H_y=0.5$ with periodic (plus imposed pressure drop of ${\rm \Delta} p=p(0,y)-p(H_x,y)=1$) and Dirichlet boundary conditions in the $x$-direction are shown in figures 5(a) and 5(b), respectively.

Figure 4.

Sketch of a slice. All fractures with centres in $[-l_f/2,H_x+l_f/2]\times [0,H_y]$ intersect with either the left or right boundary.

Figure 5.

Pressure surface plots of two realizations with ($a$) periodic and ($b$) Dirichlet boundary conditions in the $x$-direction (with an imposed pressure drop of ${\rm \Delta} p=1$). In both cases periodic boundary conditions were applied in the $y$-direction.

To solve the 1-D model equations (2.9) and (2.10), the same grid resolution in the $x$-direction was used. Next, it is explained how the model coefficients are obtained.

5.2. Model calibration

The proposed model relies on two kernel functions $\hat {g}(x,x')$ and $\hat {k}_f(x,x')$, which reflect the size, orientation and aperture distribution of the fractures, and on the effective matrix permeability $\bar {k}_m(x)$. Here, for cases with constant fracture length $l_f$, the simple ansatz

(5.1)\begin{equation} \frac{\hat{g}(x,x')}{\bar{g}} = \frac{\hat{k}_f(x,x')}{\bar{k}_f}= \begin{cases} 1, & \textrm{if }|x-x'| < l_f/2\\ 0, & \textrm{else} \end{cases} \end{equation}

is used for the kernel functions; see figure 6. This leaves us with three model constants, i.e. $\bar {k}_m$, $\bar {k}_f$ and $\bar {g}$ have to be determined in order to close the model equations (2.9) and (2.10).

Figure 6.

Normalized kernel functions $\hat {g}(x_f,x)/\bar {g}=\hat {k}_f(x_f,x)/\bar {k}_f$.

5.2.1. Effective matrix permeability

The effective matrix permeability $\bar {k}_m$ can be derived from the expression

(5.2)\begin{align} \mathbb{E}[u_m](x)&={-}\frac{k_m}{H_y}\mathbb{E}\left[\int_0^{H_y}\left((k(x,y)=k_m)\frac{\partial p}{\partial x}\right)\textrm{d}y\right]\nonumber\\ &={-}(1-\rho_f l_f a_f)k_m\frac{\partial\bar{p}_m}{\partial x} \end{align}

for the mean flow per unit area through the matrix at any location $x$. From (2.10) it becomes apparent that this flow has to match that predicted by the model, i.e.

(5.3)\begin{equation} \mathbb{E}[u_m](x)={-}\bar{k}_m\frac{\partial\bar{p}_m}{\partial x}, \end{equation}

and one obtains

(5.4)\begin{equation} \bar{k}_m=(1-\rho_f l_f a_f) k_m \end{equation}

for the effective matrix permeability. Note that the effective matrix permeability $\bar {k}_m$ is smaller than $k_m$, since some of the volume is occupied by fractures (here each fracture has a volume of $l_fa_f$). For the following studies it was assumed that $a_f$ is extremely small such that $\rho _f l_f a_f\ll 1$, and thus $\bar {k}_m$ was set to $k_m$.

5.2.2. Effective fracture conductivity

To find the effective fracture conductivity $\bar {k}_f$, flow through fractures in a very thin slice of thickness $H_x\ll l_f$ is considered; see figure 7. The average number of fractures connecting both the left and right boundaries is $\rho _f(l_f-H_x)H_y$, and therefore one can write

(5.5)\begin{equation} \mathbb{E}[u_f]\approx(l_f-H_x) \rho_fa_fk_f \frac{p_b(0)-p_b(H_x)}{H_x} \end{equation}

for the expected volumetric flow per unit area through the fractures. The model, on the other hand, predicts

(5.6)\begin{align} \mathbb{E}[u_f]&\approx\frac{p_b(0)-p_b(H_x)}{H_x}\left(\int_{H_x-l_f/2}^0\hat{k}_f(x,H_x)\, \textrm{d}x + \int_{0}^{l_f/2}\hat{k}_f(x,0)\,\textrm{d}x \right)\nonumber\\ &=(l_f-H_x)\bar{k}_f \frac{p_b(0)-p_b(H_x)}{H_x}, \end{align}

where the boundary conditions $p_b(x)=p_b(0)$ and $p_b(x)=p_b(H_x)$ were applied for all $x\le 0$ and $x\ge H_x$, respectively. Finally, one obtains the simple relation

(5.7)\begin{equation} \bar{k}_f=\rho_fa_fk_f.\end{equation}

Figure 7.

Sketch of a very thin slice. All fractures with their centres in $[H_x-l_f/2,l_f/2]\times [0,H_y]$ intersect with both the left and right boundaries.

5.2.3. Effective exchange coefficient

The effective exchange coefficient $\bar {g}$ can be determined based on the mean flow per unit area through the fractures in an infinitely large sample, which according to the model reads

(5.8)\begin{align} \mathbb{E}[u_f^\infty]&= \int_{-\infty}^0\int_0^\infty\hat{g}(x,x')\left(\bar{p}_f(x)-\bar{p}_m(x')\right) \textrm{d}x'\,\textrm{d}x\nonumber\\ &\quad +\int_{0}^\infty\int_{-\infty}^0\hat{g}(x,x')\left(\bar{p}_m(x')-\bar{p}_f(x)\right) \textrm{d}x'\,\textrm{d}x. \end{align}

Note that due to the absence of Dirichlet boundaries the term with $\hat {k}_f$ is zero. With the ansatz (5.1) for the kernel function (depicted in figure 6) this expression simplifies to

(5.9)\begin{equation} \mathbb{E}[u_f^\infty]= 2\bar{g}\int_{{-}l_f/2}^0\int_0^{x+l_f/2}\left(\bar{p}_f(x)-\bar{p}_m(x')\right) \textrm{d}x'\,\textrm{d}x \end{equation}

and with $\partial \bar {p}_m/\partial x=\nabla _x^\infty p$ one obtains

(5.10)\begin{align} \mathbb{E}[u_f^\infty]&={-}2\bar{g}\nabla_x^\infty p \int_{{-}l_f/2}^0\int_0^{x+l_f/2}\left(x'-x\right)\textrm{d}x'\,\textrm{d}x\nonumber\\ &={-}\bar{g}\frac{l_f^3\nabla_x^\infty p}{12}. \end{align}

Thus, the effective exchange coefficient can be calculated as

(5.11)\begin{equation} \bar{g}=\frac{12}{l_f^3}\frac{|\mathbb{E}[u_f^\infty]|}{|\nabla_x^\infty p|} .\end{equation}

To obtain $\mathbb {E}[u_f^\infty ]$, a Monte Carlo study with periodic boundary conditions plus imposed mean pressure gradient of $\nabla _x^\infty p=-H_x^{-1}$ in the $x$-direction was performed. Domain size and grid resolution were chosen as described in § 5.1, i.e. $H_x=1.5$, $H_y=0.5$ and $h_x=h_y=0.01$. Numerical results of 1000 independent realizations were sampled, which resulted in $|\mathbb {E}[u_f^\infty ]|/|\nabla _x^\infty p|\approx 4.285$ with a statistical standard error of $0.1\,\%$.

As a result, the model coefficients used for the following comparative validation study are $\bar {k}_m=1$, $\bar {k}_f=100$ and $\bar {g}=3290.88$.

5.3. Model validation

In order to obtain reference data for the model validation, Monte Carlo studies with Dirichlet boundary conditions at $x=0$ and $x=H_x$, that is, with $p_b(0,y)=1$ and $p_b(H_x,y)=0$, and periodic boundary conditions in the $y$-direction were performed. Note that for the individual realizations of the Monte Carlo studies the permeability in the cells along the horizontal pixel lines representing fractures was set to $k_f^{numerical}=k_fa_f/h_y=200$ and in the cells representing the porous matrix it was $k_m^{numerical}=1$. However, the way the homogenized equations were derived is independent of the ratio $k_f^{numerical}/k_m^{numerical}$ and the model is not limited to high contrast or low contrast cases. Pressure solutions were computed for 5000 independent permeability fields, where each realization was generated as described in § 5.1. In all cases the agreement between model solutions and Monte Carlo is excellent. Note that the sharp steps in $\bar {p}_f$ at a distance of $l_f/2$ from the boundaries are characteristic for cases with a single fracture length, since all fractures within that distance are directly connected to the boundaries, while for all others the boundary pressure signal first has to bass through the matrix. Figure 8 shows profiles of expected matrix pressure $\bar {p}_m(x)$ (a,c,e) and expected fracture pressure $\bar {p}_f(x)$ (b,d,f) for domains with $H_x\in \{0.3,0.6,0.9\}$; symbols refer to Monte Carlo data and the solid lines to model solutions; the dotted lines will be addressed later. Note that $\bar {p}_f(x)$ refers to the expected pressure of fractures with their centres at location $x$. Figure 9(a) shows the normalized flows $\mathbb {E}[u](H_x)H_x$ obtained from Monte Carlo and the non-local model as functions of the slice thickness $H_x$. As expected, the effective conductance decreases with larger $H_x$; a similar behaviour was observed at a smaller scale for pore networks (Meyer & Gomolinski Reference Meyer and Gomolinski2019). Note, however, that they do not distinguish between different media, but applied the same framework as Delgoshaie et al. (Reference Delgoshaie, Meyer, Jenny and Tchelepi2015) and demonstrate how the parameters can be extracted from pore network data. For each $H_x\in \{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.5\}$, 1000 independent solutions were sampled; the statistical standard error of the normalized flows obtained from Monte Carlo is shown in figure 10. It can be observed that for very small and large $H_x$-values the model predictions compare very well with the Monte Carlo data. The discrepancies for $H_x\approx l_f$ can be attributed to the use of a constant $\bar {g}$, the ansatz with piecewise constant kernel functions and of course the intrinsic simplifications made for the model derivation. Next, the effect of the former is investigated. It is argued that fluid exchange between fractures and matrix near the Dirichlet boundaries is underestimated (in comparison to the domain centre), which suggests using higher $\bar {g}$ values in the vicinity of such domain boundaries. Similar non-local effects were observed in pore network studies Meyer & Gomolinski (Reference Meyer and Gomolinski2019). The approach followed here to correct boundary effects not accounted for by the second right-hand-side term in (2.8) is inspired by elliptic relaxation models used to describe turbulent flows. There, the goal is a similar one, i.e. to distinguish between the flow dynamics near walls and in free shear flow regions; mainly due to wall blocking effects. Elliptic relaxation models can account for such non-local wall effects in a very elegant way. In the case considered here, the constant $\bar {g}$ value is replaced by $\tilde {g}(x)$, which is the solution of the elliptic differential equation

(5.12)\begin{equation} \left(1-\left(\frac{l_f}{2}\right)^2\frac{\partial^2}{\partial x\partial x}\right)\tilde{g}=\bar{g}, \end{equation}

with the Dirichlet boundary condition $\tilde {g}(0)=\tilde {g}(H_x)=\tilde {g}^{boundary}$. In the statistically one-dimensional case, its analytic solution is

(5.13)\begin{equation} \tilde{g}(x)= \bar{g} + \frac {(\tilde{g}^{boundary}-\bar{g})\left({\exp}({2x/l_f})+\exp({2(H_x-x)/l_f})\right)} {1+\exp({2H_x/l_f})}, \end{equation}

but in general, e.g. for most 2-D and 3-D problems, (5.12) has to be solved numerically. Note that (5.12) ensures that $\tilde {g}(x)$ relaxes towards $\bar {g}$ away from the boundaries, whereas the spatial correlation length scale is of the order of $l_f/2$, which seems to be an obvious choice. The elliptic equation (5.12) also reflects the elliptic nature of the flow problem (see (2.5)). For the cases investigated so far, it was empirically found that the best match with Monte Carlo data is obtained with $\tilde {g}^{boundary}=1.25\bar {g}$. Figure 9(b) shows the obtained normalized flow rates as a function of $H_x$, and it can be observed that the match with the Monte Carlo reference data is excellent. One can also see that the corresponding pressure profiles shown in figure 8 (dashed lines) are in slightly better agreement with the Monte Carlo reference than those obtained without elliptic relaxation. Note, however, that elliptic relaxation only leads to a noticeable improvement if $H_x\approx l_f$; otherwise, the agreement is already very good without elliptic relaxation.

Figure 8.

(a,c,e) Profiles of expected matrix pressure $\bar {p}_m(x)$; (b,d,f) expected fracture pressure $\bar {p}_f(x)$ (representing the pressures of fractures with their centres at location $x$) for $H_x\in \{0.3,0.6,0.6\}$. Lines represent model predictions and symbols Monte Carlo data.

Figure 9.

Normalized flow $\mathbb {E}[u]H_x$ computed with the non-local model for different domain sizes in comparison with the Monte Carlo data. Left without and right with elliptic relaxation of $\tilde {g}(x)$.

Figure 10.

Statistical standard error of the normalized flow $\mathbb {E}[u]H_x$ obtained from Monte Carlo with 1000 realizations for each $H_x$-value.

Regarding the cost of solving the homogenized equations, one has to distinguish between preprocessing and the actual simulations; the former is necessary to obtain closed expressions for $\hat {g}$, $\hat {k}_f$ and $\bar {k}_m$, and the latter to obtain $\bar {p}_m$ and $\bar {p}_f$. Compared to Monte Carlo solving the homogenized integral equations with a known parameter set is very cheap, that is, only one realization has to be computed and only mean variations need to be spatially resolved. The only extra cost results from the discretization of the integral terms, which leads to a broader stencil and a linear system with a fuller matrix. The preprocessing step typically required Monte Carlo with a stationary (statistically homogeneous) setting, but once determined the same kernels can be used for any non-stationary setting with the same statistical fracture distribution. For the cases discussed in the paper the kernel shape (for both $\hat {g}$ and $\hat {k}_f$) is very simple; as shown in figure 6. This shape will change, however, if more complicated fracture length, aperture and orientation distributions are to be considered, which is a subject of future research. Here, the magnitude of $\hat {k}_f$, i.e. $\bar {k}_f$ is analytically determined according to (5.7). The magnitude of $\hat {g}$, i.e. $\bar {g}$, is determined based on (5.11), which requires us to empirically determine $\mathbb {E}[u_f^\infty ]$. Therefore, a Monte Carlo study (consisting of $1000$ realizations) with an imposed mean pressure gradient in the $x$-direction of $\nabla _x^\infty p=-H_x^{-1}$ was performed.