The first integral vanishes because $\nabla \cdot \mathbf{v}' = 0$, and the second integral can be identified as the volume average flux $\langle \mathbf{v}' \rangle_{VAF}$ (see Table 2). Therefore, the volume average of equation (Equation (99)) is equal to

$$-\frac{1}{\Omega} \int_{\Omega} \nabla h \cdot (\mathbf{k} \nabla h') d\omega = -\mathbf{J} \cdot \mathbf{v}'_{VAF} = -\mathbf{J} \cdot (\mathbf{K}_{VAF} \mathbf{J}') \quad (151)$$

The previous result is obtained using the definition of $\mathbf{K}_{VAF}$, and the fact the volume average gradient of $h'$ is equal to $-\mathbf{J}'$, in addition to Darcy's law.

Finally, exactly the same averaging and derivation can be made on the right-hand side of Equation (147), and therefore we obtain the following identity, implying that $\mathbf{K}_{VAF}$ is symmetric.

$$\mathbf{J} \cdot (\mathbf{K}_{VAF} \mathbf{J}') = \mathbf{J}' \cdot (\mathbf{K}_{VAF} \mathbf{J}) \quad (152)$$

### 7.4. VAF + GRD Is Positive Definite

$\mathbf{K}_{VAF}$ is symmetric when it is computed with GRD, but it is also positive and definite. This is a direct consequence of the fact that under the GRD boundary conditions, the VAF criteria ensures that the Total Dissipated Power is identical in the upscaled medium and in the heterogeneous medium. Indeed, one can show that $\mathbf{K}_{VAF}$ is positive and definite if for any vector $\mathbf{J}$ different from 0 the following quadratic form is always positive:

$$\mathbf{J} \cdot \mathbf{K}_{VAF} \mathbf{J} > 0 \quad (153)$$

This is always true because for whatever $\mathbf{J}$, the expression above is proportional to the Total Dissipated Power in the medium if we apply linearly varying head boundary conditions (GRD conditions). Since the small-scale hydraulic conductivity tensor $\mathbf{k}$ is assumed positive definite, the total dissipated power on the heterogeneous medium (left hand side of equation below) is strictly positive when the gradient imposed to the block is different from zero.

$$P = -\rho_g \mathbf{J} \cdot \langle \mathbf{v} \rangle_{VAF} = \rho_g \mathbf{J} \cdot (\mathbf{K}_{VAF} \mathbf{J}) \quad (154)$$

The right-hand side of the above equation is the Total Dissipated Power on the upscaled medium, which is equal to the TDP on the homogeneous medium. The fluid density and the acceleration of gravity are strictly positive; the right-hand side of the previous equation is then strictly positive, and $\mathbf{K}_{VAF}$ is positive definite.

### 7.5. Symmetry of VAF / TDP & PRD

Concerning the symmetry property of the equivalent conductivity under VAF averaging for the TDP criterion and the permeametric PRD criterion, Boe [96] and Farmer [35] used Equation (72) to define the equivalent conductivity tensor and prove that the resulting tensor is symmetric. The result is in agreement with the numerical experiments of Durlofsky [37].

## 8. Discussion and Conclusions

It was seen in the preceding review and analyses that several other works in the literature have compared upscaling techniques. For instance, a large comparison exercise of upscaling techniques has been made by the Society of Petroleum Engineers. The problem included the simulation of multiphase flow in a large heterogeneous oil reservoir after applying various upscaling techniques. The results were reported by Christie and Blunt [108]. One of the findings of this exercise was that the equivalent conductivity obtained with no flow boundary conditions were providing "better" upscaled solutions than the results obtained with linearly varying heads (or gradient conditions) for their specific case study.

Here in this paper, we take a different perspective by conducting a more theoretical analysis of various approaches that have been used in the past to define the equivalent permeability tensor. We started by setting the problem in a very open manner and trying to review a broad range of definitions and practical ways to compute the equivalent permeability.
ability on finite domains or "blocks". We then analyzed and illustrated the various methods with simple numerical experiments and with analytical proofs concerning their properties.

More precisely, the definitions of the equivalent conductivity that we considered involved different boundary conditions such as permeametric conditions, linear head conditions (or fixed gradient condition), and different equivalence criteria such as Volume Average Flux, Net Surface Flux, Vectorial Surface Flux, and Total Dissipated Power). We show that not only does the value of the computed equivalent permeability depends on the type of averaging technique and boundary conditions, but also the properties of the resulting equivalent conductivity can be different. These results were summarized in Table 4.

For a medium stratified along the axes of a rectangular block, permeameter type boundary conditions (PRM) allow to obtain the theoretical equivalent conductivity (which coincides in that case with the theoretical infinite domain "effective" conductivity). In contrast, for the same case of stratification parallel to the block's axes, linearly varying head conditions (gradient conditions "GRD") overestimate the smallest component of the equivalent conductivity tensor. Indeed, the permeameter imposes confining conditions that imply a no flow boundary perpendicularly to the strata, whereas linear head conditions generate small scale lateral circulation cells within the high conductivity strata, which increase flux perpendicular to strata (see Figure 16) leading to an overestimation of total flux. On the other hand, when the strata are inclined with respect to the block's axes, linearly varying head conditions provide a better estimate, closer to the expected result than with permeameter type conditions.

Among the various criteria that can be used, the Volume Average Flux (VAF) criterion is particularly interesting because of the properties of the resulting equivalent permeability tensor:

* $K_{ij}$ is always a symmetric tensor when combined with GRD (prescribed gradient), FLX (prescribed flux), or PRD (periodic) boundary conditions;

* The VAF method using GRD, FLX, or PRD boundary conditions, ensures that the Total Dissipated Power (TDP) is identical in the heterogeneous medium and in the equivalent homogeneous medium.

The numerical experiments show that non-symmetry of the equivalent $K_{ij}$ tensor occurs more frequently with other equivalence criteria (NSF & VSF).

Negative values of some of the principal permeability values $K_{ii}^*$ of the Kij tensor can occur with the VSF equivalence criterion. Recall that the principal values $K_{ii}^*$ are the eigenvalues of Kij, and these are independent of the chosen coordinate system. With the VSF criterion, it is found that the resulting equivalent Kij tensor can be indefinite (e.g., with one negative principal value and the others positive), although it cannot be negative-definite (which is a physically sound result).

In addition to continuous and composite porous media, we have also presented partial results on equivalent permeability for some simple 2D fracture networks (Sections 5.1 and 6.2). The algebraic properties of equivalent permeability of irregular fracture networks, represented as networks of links in 2D space, have also been analyzed by Renard and Ababou [62], Ababou and Renard [61], providing partial results for at least some of the criteria and conditions reviewed in the present paper. A more complete analysis is under way [97]. The case of networks of planar fractures in 3D space is another matter, which involves a richer class of "links" and "intersections" compared to 2D networks. We refer the reader to Noetinger and Jarrige [109] for semi-analytical/numerical approaches to 3D fracture network flow, to Rajeh et al. [64] for a semi-empirical approach to the equivalent permeability tensor of 3D fractured porous media with permeable porous matrix, and to Cañamón et al. [63], and references therein, for a generalized graph-based approach to percolation phenomena in 3D networks of planar fractures.

Recalling the possibility of using flux-based or velocity-based flow equations (Sections 2.2.2 and 3.1, Appendix A.2), we show that the averaged version of this flow equation may not yield the same equivalent permeability than the head-based flow equa-
tions. For instance, the velocity version of Darcy's law, Equation (7), may lead, upon averaging, more naturally, to an equivalent resistivity rather than an equivalent permeability. The challenge may then be to match, or at least to compare, the equivalent permeability from head-based equations to equivalent resistivity from velocity-based equations. A "dual" permeability-resistivity matching procedure (summarized in Section 4.3) was presented by Fadili and Ababou [75] for upscaling in the context of randomly heterogeneous media (effective properties in an infinite domain). A more general matching procedure, which we leave for future work, could be attempted for arbitrarily heterogeneous porous and/or fractured finite blocks, in combination with the various averaging operations and boundary conditions that were tested in the present work.

In closing, to sum up, we have presented in this paper a comprehensive review as well as theoretical and numerical analyses on upscaling methods and criteria, aimed at calculating equivalent block permeability tensors of heterogeneous media, and also, at evaluating the properties of the proposed upscaling criteria.

**Author Contributions:** Writing—original draft preparation, P.R., R.A.; writing—review and editing, P.R., R.A. Both authors contributed equally to the article and agreed to its publication. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A.

### Appendix A.1. Tensorial Kij, Directional K, and Anisotropic Kij Ellipse

This Appendix presents a more detailed algebraic complement to Section 2.4 in the text, concerning second rank tensors such as the conductivity tensor $K_{ij}$, and several related concepts, including: principal axes, principal values, directional "flux" and "gradient" conductivities, anisotropic conductivity ellipses.

#### Appendix A.1.1. Vectors, Second Rank Tensors, and Tensorial Conductivity

Vectors, Tensors, and Transformation Rules

A vector is a first rank tensor. A second rank tensor $T$ or $T_{ij}$ ($i = 1, \dots, N; j = 1, \dots, N$) can be represented by a matrix, and it can be defined as follows, by considering that the tensor $T$ is a linear application from $\mathbb{R}^N$ to $\mathbb{R}^N$:

$$ T: \mathbb{R}^N \rightarrow \mathbb{R}^N; \mathbf{u} \rightarrow \mathbf{v} = T(\mathbf{u}) = \mathbf{Tu} \quad (A1) $$

or: $u_i \rightarrow v_i = T_{ij}u_j$ (with implicit summation on repeated indices)

In 2D, with $N = 2$, the application $T$ is also represented by a square $2 \times 2$ matrix:

$$ T \text{ as an application } \begin{bmatrix} v_1 \\ v_2 \end{bmatrix} = \begin{bmatrix} T_{11} & T_{12} \\ T_{21} & T_{22} \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \end{bmatrix}; \quad T \text{ as a square matrix } T = \begin{bmatrix} T_{11} & T_{12} \\ T_{21} & T_{22} \end{bmatrix} \quad (A2) $$

The coefficients ($T_{ij}$) depend on the system of coordinates $(\mathbf{e}_1, \mathbf{e}_2)$, and the $T_{ij}$'s change if the system is changed by a rotation $\mathbf{R}: (\mathbf{e}_1, \mathbf{e}_2) \rightarrow (\mathbf{e}_1', \mathbf{e}_2')$. Conversely, its transpose $\mathbf{P} = \mathbf{R}^T$ transforms the new system into the initial system, and is called the passage matrix $\mathbf{P}: (\mathbf{e}_1', \mathbf{e}_2') \rightarrow (\mathbf{e}_1, \mathbf{e}_2)$. In 2D, we have:

$$ \mathbf{R} = \begin{bmatrix} \cos(\theta) & +\sin(\theta) \\ -\sin(\theta) & \cos(\theta) \end{bmatrix} \quad \text{and} \quad \mathbf{P} = \mathbf{R}^T = \begin{bmatrix} \cos(\theta) & -\sin(\theta) \\ +\sin(\theta) & \cos(\theta) \end{bmatrix} \quad (A3) $$

The rotation matrix and the passage matrix are orthogonal: $\mathbf{P}\mathbf{P}^T = \mathbf{P}^T\mathbf{P} = \mathbf{I}$ (identity matrix).
In 3D, the passage matrix P becomes:

$$ \mathbf{P} = \mathbf{R}^{\mathrm{T}} = \begin{bmatrix} \cos(\theta) \cos(\varphi) & -\sin(\theta) & -\cos(\theta) \sin(\varphi) \\ +\sin(\theta) \cos(\varphi) & \cos(\theta) & -\sin(\theta) \sin(\varphi) \\ \sin(\varphi) & 0 & \cos(\varphi) \end{bmatrix} \quad (\mathrm{A}4) $$

The transformation rule for a vector **u** under a change of coordinate system, that is, under a rotation of the basis vectors (**e**₁, **e**₂) → (**e**₁', **e**₂'), is as follows, in terms of the passage matrix **P**:

$$ \mathbf{u} \rightarrow \mathbf{u}' = \mathbf{P}^{\mathrm{T}} \mathbf{u} \mathbf{u}' \rightarrow \mathbf{u} = \mathbf{P} \mathbf{u}' \quad (\mathrm{A}5) $$

Note: vectors **u** and **u'** can represent the position vectors (**x** and **x'**) in the old and new system, respectively.

The quantity **u** is a "true" vector only if it is transformed as shown above under a rotation of the coordinate system. Otherwise, it is not a true vector. Example: if we have **u** = [1 1] in the reference system and **u'** = [1 0] in the new system rotated by angle θ = π, then **u** is not a true vector because it does not obey the transformation rule (it is, rather, a pseudo-vector).

Finally, the transformation rule for a 2nd rank tensor **T** under a rotation of the basis vectors (**e**₁, **e**₂) → (**e**₁', **e**₂') is as follows, in terms of the passage matrix **P**:

$$ \mathbf{T} \rightarrow \mathbf{T}' = \mathbf{P}^{\mathrm{T}} \mathbf{T} \mathbf{P}; \quad \mathbf{T}' \rightarrow \mathbf{T} = \mathbf{P} \mathbf{T}' \mathbf{P}^{\mathrm{T}} \quad (\mathrm{A}6) $$

$$ T_{ij} \rightarrow T'_{kl} = P_{ik}T_{ij}P_{jl}; \quad T'_{kl} \rightarrow T_{ij} = P_{ik}T'_{kl}P_{jl} \quad (\text{summation on repeated indices}) $$

## Application to Darcy's Law with Anisotropic $K_{ij}$ Tensor

Darcy's law is $\mathbf{v} = -\mathbf{K} \nabla h$, where the Darcy velocity vector $\mathbf{v}$, the head gradient vector $\nabla h$, and the conductivity tensor $\mathbf{K}$, are all expressed in the reference system ($\mathbf{e}_i$). The notation $\mathbf{J} = -\nabla h$ is also used to designate the hydraulic gradient and will be used elsewhere.

We now derive the expression of Darcy's law in a new transformed system ($\mathbf{e}_i'$) as follows. We start with Darcy's law in the reference system ($\mathbf{e}_i$) and we insert the vector transformation rules $\mathbf{v} = \mathbf{P} \mathbf{v}'$ and $\nabla h = \mathbf{P} \nabla h'$, and then we also insert the tensor transformation rule $\mathbf{K} = \mathbf{P} \mathbf{K}' \mathbf{P}^{\mathrm{T}}$ for the conductivity. This yields:

$$ \mathbf{v}' = - (\mathbf{P} \mathbf{K}' \mathbf{P}^{\mathrm{T}}) \mathbf{P} \nabla h' = -\mathbf{P} \mathbf{K}' \nabla h' \Rightarrow \mathbf{v}' = -\mathbf{K}' \nabla h' \quad (\mathrm{A}7) $$

which is just Darcy's law in the new system ($\mathbf{e}_i'$), with the transformed conductivity $\mathbf{K}'$ as expected. This confirms, or proves, that the permeability tensor is indeed transformed by the rule $\mathbf{K} = \mathbf{P} \mathbf{K}' \mathbf{P}^{\mathrm{T}}$. We also have conversely $\mathbf{K}' = \mathbf{P}^{\mathrm{T}} \mathbf{K} \mathbf{P}$, using the fact that $\mathbf{P} \mathbf{P}^{\mathrm{T}} = \mathbf{P}^{\mathrm{T}} \mathbf{P} = \mathbf{I}$.

## Principal System, Diagonalization

These transformation rules can be used in particular for finding the principal vectors ($\mathbf{e}_i^*$) and the principal permeability values ($K_{ii}^*$ without summation). In the new principal system ($\mathbf{e}_i^*$), the conductivity tensor $K_{ij}$ is expressed by a diagonal matrix diag($K_{11}^*, K_{22}^*, K_{33}^*$) with real non-negative values. Briefly, the principal vectors (eigenvectors) and the principal values (eigenvalues) of the conductivity tensor $K_{ij}$ are obtained by solving the equations:

$$ \det(\mathbf{K} - \lambda \mathbf{I}) = 0 \Rightarrow \lambda = \begin{vmatrix} \lambda_1 \\ \lambda_2 \\ \lambda_3 \end{vmatrix} \quad (\mathrm{A}8) \\ \text{then: } \begin{vmatrix} \mathbf{K} \mathbf{e}_1^* = \lambda_1 \mathbf{e}_1^* \\ \mathbf{K} \mathbf{e}_2^* = \lambda_2 \mathbf{e}_2^* \\ \mathbf{K} \mathbf{e}_3^* = \lambda_3 \mathbf{e}_3^* \end{vmatrix} = [\mathbf{e}_1^*, \mathbf{e}_2^*, \mathbf{e}_3^*] $$
The first equation has three positive real roots ($\lambda_1, \lambda_2, \lambda_3$) if $\mathbf{K}$ is symmetric positive definite; the second equation consists in solving for the three principal vectors $[\mathbf{e}_1^*, \mathbf{e}_2^*, \mathbf{e}_3^*]$ once the eigenvalues ($\lambda_1, \lambda_2, \lambda_3$) are known. The principal basis vectors ($\mathbf{e}_i^*$), or eigenvectors, are the column vectors of the passage matrix $\mathbf{P}^*: (\mathbf{e}_i^*) \to (\mathbf{e}_i)$ that takes the new principal basis $\mathbf{e}_i^*$ into the initial basis $\mathbf{e}_i$ (such that $\mathbf{e}_i = \mathbf{P}^* \mathbf{e}_i^*$). In the new principal basis, the conductivity tensor $K_{ij}$ is represented by the diagonal matrix:

$$ \mathbf{K}^* = \begin{bmatrix} K_{11}^* & 0 & 0 \\ 0 & K_{22}^* & 0 \\ 0 & 0 & K_{33}^* \end{bmatrix} \quad \text{where } K_{11}^* = \lambda_1, \ K_{22}^* = \lambda_2, \ K_{33}^* = \lambda_3 \quad (A9) $$

The conductivity matrix $\mathbf{K}$ in a given coordinate system can therefore always be expressed in terms of the diagonalized matrix $\mathbf{K}^*$ using the following relation, where the $\mathbf{e}_i^*$'s are column vectors:

$$ \mathbf{K} = \mathbf{P}^* \mathbf{K}^* \mathbf{P}^* \mathbf{T} = \mathbf{P}^* \begin{bmatrix} K_{11}^* & 0 & 0 \\ 0 & K_{22}^* & 0 \\ 0 & 0 & K_{33}^* \end{bmatrix} \mathbf{P}^* = [\mathbf{e}_1^*, \mathbf{e}_2^*, \mathbf{e}_3^*] \begin{bmatrix} K_{11}^* & 0 & 0 \\ 0 & K_{22}^* & 0 \\ 0 & 0 & K_{33}^* \end{bmatrix} [\mathbf{e}_1^*, \mathbf{e}_2^*, \mathbf{e}_3^*]^T \quad (A10) $$

It is advisable to sort the eigenvalues (principal conductivities) in decreasing order, i.e., $K_{11}^* \ge K_{22}^* \ge K_{33}^*$.

Thus, in 2D, let “$\gamma$” designate the angle between the first principal axis $\mathbf{e}_1^*$ (associated with the first and largest principal conductivity $K_{11}^*$), and the first axis $\mathbf{e}_1$ of the reference system. This angle $\gamma = (\mathbf{e}_1^*, \mathbf{e}_1)$ can be determined to within an additive constant $\pi$ from the relation $\tan(\gamma) = (K_{11}^* - K_{11})/K_{12}$. If the eigenvalues are not ordered, here is another relation that determines $\gamma$ independently of the ordering of principal values, but only within an additive constant $\pi/2$: $\tan(2\gamma) = K_{12}/(K_{11} - K_{22})$. The 2D principal conductivity values are given by:

$$ K_{11}^* = \frac{K_{11} + K_{22}}{2} + \left\{ \left( \frac{K_{11} - K_{22}}{2} \right)^2 + (K_{12})^2 \right\}^{1/2} \quad (A11) $$

$$ K_{22}^* = \frac{K_{11} + K_{22}}{2} - \left\{ \left( \frac{K_{11} - K_{22}}{2} \right)^2 + (K_{12})^2 \right\}^{1/2} $$

In particular, we retrieve the identity $K_{11}^* + K_{22}^* = K_{11} + K_{22}$, as it should, since the trace of a tensor $K_{ij}$ is invariant under a change of the coordinate system.

Note: the above description concerns the diagonalization of a *symmetric positive definite* $K_{ij}$ tensor (For a non-symmetric tensor, there are two distinct principal bases and two distinct sets of principal values.).

## Appendix A.1.2. Directional Conductivity, and Conductivity Ellipse

If $K_{ij}$ is a second rank tensor, symmetric and positive definite (or, more generally, this applies also to non-symmetric $K_{ij}$ if its symmetric part is definite positive), then it can also be described equivalently in the form of a directional conductivity parallel to flux ($K_{FLUX}$), or parallel to hydraulic gradient ($K_{GRAD}$). They are functions of the direction vector of Darcy velocity $\mathbf{v}$ or of the gradient $\mathbf{J} = -\nabla h$, respectively. The direction angles can be used to indicate the direction; polar angle $\theta$ in 2D; spherical angles ($\theta, \phi$) in 3D. It can then be shown that the polar plots of these directional conductivities (or their square roots) describe an ellipse in 2D, or an ellipsoid in 3D. The conductivity ellipse contains all the information about the $K_{ij}$ tensor: the principal radii are related to the principal values of $K_{ij}$, and the principal axes of the ellipse coincide with the principal vectors of $K_{ij}$.
Directional Conductivity

Let us start again with the tensorial form of Darcy's law, using the notation $\mathbf{J} = -\nabla h$:

$$ \mathbf{v} = -\mathbf{K} \nabla h = \mathbf{K} \mathbf{J} \qquad (A12) $$

or $v_i = K_{ij} J_j$ (with summation on repeated indices)

Let $\mathbf{n}_F$ be the unit “flux” vector (aligned with Darcy velocity); then, by definition: $\mathbf{n}_F = \mathbf{v}/\|\mathbf{v}\|$. Similarly, let $\mathbf{n}_G$ be the unit “gradient” vector (aligned with $\mathbf{J} = -\nabla h$); then by definition: $\mathbf{n}_G = \mathbf{J}/\|\mathbf{J}\|$. The directional “flux” conductivity ($K_{FLUX}$) is a scalar quantity obtained by “projecting” Darcy’s law on the direction of the “flux”, as follows:

$$ K_{FLUX} \text{ is such that: } \|\mathbf{v}\| = K_{FLUX} (\mathbf{J} \cdot \mathbf{n}_F) \Rightarrow K_{FLUX} = \|\mathbf{v}\|^2 / (\mathbf{J} \cdot \mathbf{v}) \qquad (A13) $$

Similarly, the directional “gradient” conductivity ($K_{GRAD}$) is a scalar quantity obtained by “projecting” Darcy’s law on the direction of the “gradient”, as follows:

$$ K_{GRAD} \text{ is such that: } (\mathbf{v} \cdot \mathbf{n}_G) = K_{GRAD} \|\mathbf{J}\| \Rightarrow K_{GRAD} = (\mathbf{v} \cdot \mathbf{J}) / \|\mathbf{J}\|^2 \qquad (A14) $$

Anisotropic Conductivity Ellipse

By re-inserting Darcy’s law with tensorial conductivity, and noting that $\|\mathbf{v}\|^2 = \mathbf{v}\mathbf{v}^T$ and $\|\mathbf{J}\|^2 = \mathbf{J}\mathbf{J}^T$, it can be shown that [40,77]:

$$ K_{FLUX} = (\mathbf{n}_F^T \mathbf{K}^{-1} \mathbf{n}_F)^{-1} K_{GRAD} = (\mathbf{n}_G^T \mathbf{K}^{-1} \mathbf{n}_G)^{-1} \quad \text{and} \quad K_{GRAD} = \mathbf{n}_G^T \mathbf{K} \mathbf{n}_G \qquad (A15) $$

In 2D, let $\theta_F$ and $\theta_G$ designate the polar angles of the “flux” (Darcy velocity) and of the “gradient” (hydraulic gradient). Then, the direction vector $\mathbf{n}$ is of the form $\mathbf{n} = [\cos(\theta), \sin(\theta)]$, and the directional conductivities given just above are each a function of the polar angle $(\theta): K_{FLUX}(\theta_F)$ and $K_{GRAD}(\theta_G)$.

Now, it can be shown by inserting the previous diagonalization $\mathbf{K} = \mathbf{P}^* \mathbf{K}^* \mathbf{P}^* \mathbf{T}$ in these expressions that each of the scalar quantities, $\sqrt{K_{FLUX}(\theta_F)}$ and $1/\sqrt{K_{GRAD}(\theta_G)}$, are the polar radii describing an ellipse whose principal axes are the eigenvectors of $\mathbf{K}$, and whose principal radii are $\sqrt{K_{ii}^*}$ (for the directional conductivity // flux ellipse) and $1/\sqrt{K_{ii}^*}$ (for the directional gradient // ellipse). Similarly, in 3D, $\sqrt{K_{FLUX}(\theta_F, \phi_F)}$ and $1/\sqrt{K_{GRAD}(\theta_G, \phi_G)}$ are the polar radii describing an ellipsoid for the “flux” and “gradient” versions, respectively.

Many authors have applied these concepts to analyze the (upscaled) hydraulic conductivity of heterogeneous or fractured media. For instance, Long et al. (1982) designed 2D numerical flow experiments under imposed hydraulic gradient $\mathbf{J}_0$ in a synthetic fractured rock; they computed the directional $K_{GRAD}$ along the gradient direction, and then rotated $\mathbf{J}_0$ to obtain $K_{GRAD}$ for different angular directions, plotted the resulting $1/\sqrt{K_{GRAD}(\theta_G)}$, and fitted the resulting plot to an ellipse. Applications of the ellipse representation to equivalent block conductivity anisotropy can also be found in the main text: Figures 15 and 16.

Finally, note that the symmetric part of the $K_{ij}$ tensor, $\mathbf{K}_S = (\mathbf{K} + \mathbf{K}^T)/2$, can be analyzed in the same way as above (i.e., in terms of directional conductivities and ellipses). This is because the quadratic form $\mathbf{x}^T \mathbf{K} \mathbf{x}$ of $\mathbf{K}_S$ and the non-symmetric $\mathbf{K}$ are the same. However, the anti-symmetric part, $\mathbf{K}_A = (\mathbf{K} - \mathbf{K}^T)/2$, has a zero quadratic form, and the physical role of $\mathbf{K}_A$ cannot be analyzed in terms of directional conductivities or in terms of conductivity ellipses. The complex case of non-symmetric conductivity anisotropy is discussed in the main text around Figure 2 in the context of equivalent conductivity models.
## Appendix A.2. Derivation of Velocity-Based PDE Governing Darcy Flow

### Appendix A.2.1. Introduction to Flux Based or Velocity Based Flow Equation

In this Appendix, we provide a brief derivation of the flux-based PDE governing the Darcy velocity vector (rather than the pressure or hydraulic head) in a 3D heterogeneous porous medium, after the work of Ababou [68] (Section 4.3, pp. 211–217), which is also presented in Zijl and Nawalany [69] (Section 3.5 “Directly calculated velocity components”, pp. 55–63).

First, this development is motivated by the fact that the spatial structure of the Darcy velocity vector is important, since it has a direct influence on the spatial distribution, variability, and dispersion of advected tracers and contaminants, and also, since it plays a major role in the determination of effective permeability (the main topic of this paper). Note that pore water velocity (which determines solute advection) is directly related to the Darcy velocity or flux vector, since the pore-velocity may be defined locally as flux divided by porosity.

NB: to be complete, note also that we use the term “flux” for “areal flux density” denoted $Q$, the latter being also identical to the so-called “Darcy velocity” denoted $V$ (thus we use the terms “flux” or “Darcy velocity” interchangeably, and the notations $V$ and $Q$ stand for the same vector quantity).

The flux-based equation was used in the literature to develop perturbation solutions for the flux vector in randomly heterogeneous media, as an alternative to the indirect determination of flux obtained by inserting head perturbation solutions in an approximate (perturbed) form of Darcy’s law. The more direct approach based on flux-governing equations has led to first order (and possibly higher order) stochastic flux solutions. This direct method of stochastic flux analysis was developed earlier in a Ph.D. thesis (Ababou 1988, Chap. 4, Section 4.3), which is also available as an MIT report 1988 by Ababou, McLaughlin, Gelhar. In addition, this flux approach has also been later presented and discussed in detail in Zijl and Nawalany’s [69] book (Chap.3) directly based on the aforementioned work. Some results from this direct stochastic flux equation were presented in [11].

Here, we present the general flux-based equation, or velocity-based equation (in the sense of Darcy velocity) for heterogeneous porous media (random or not), without presenting the stochastic flux analysis dedicated specifically to random auto-correlated continua.

### Appendix A.2.2. Development of Velocity Vector Equations

At steady state, the flux or Darcy velocity vector $v$ is governed by mass conservation $div(v) = 0$ and Darcy’s law $v = -K \text{grad}(h) = K J$ where $J$ is the “hydraulic gradient” vector defined as $J = -\text{grad}(h)$. Note that, by Darcy’s law, the hydraulic gradient obeys $J = v/K$ where $K = K(x)$ for a heterogeneous porous medium (here we consider locally isotropic heterogeneity, so that $K(x)$ is a scalar field). Now, mass conservation implies that $v$ is divergence free (solenoidal), and Darcy’s law implies that the hydraulic gradient $J = v/K$ is curl free (irrotational). Let us now define the vorticity vector $\Omega$ as the curl of velocity, that is:

$$ \Omega = \nabla \times v \qquad (A16) $$

Inserting $v = K J$ and using the fact $\nabla \times J = 0$ (since the curl of a gradient is always null) yields a relation between vorticity, log-conductivity, and flux, as follows:

$$ \Omega = \nabla LnK \times v = \nabla F \times v \qquad (A17) $$

where we denote $F(x)$ the log-conductivity field: $F(x) = LnK(x)$.

On the other hand, since $v$ is divergence free, the gradient of its divergence is the null vector. Developing the gradient of divergence by the standard rules of vector field operators (e.g., Gradshtein and Rhyzik, 2007, [110], Section 10.31 therein): Properties of the Operator $\nabla$) yields finally:

$$ \nabla^2 v + \nabla \times \Omega = 0 \qquad (A18) $$
Note that this is a vector equation. The operator $\nabla^2$ stands for the vector Laplacian, i.e., the Laplacian operator applied to each component of vector $\mathbf{v}$.

Finally, substituting Equation (A17) into Equation (A18) yields the vectorial velocity-based flow equation:

$$ \nabla^2 \mathbf{v} + \nabla \times (\nabla F \times \mathbf{v}) = \mathbf{0} \qquad (A19) $$

Following Ababou ([68], Section 4.3 therein), this velocity-equation (also called “flux equation”) can be further developed by decomposing the term $\nabla \times (\nabla F \times \mathbf{v})$ using the standard rules of vector field operators [110]. However, before proceeding with this decomposition, the vector equation can be developed in a more instructive way by introducing physically meaningful quantities, analogous to those used in elasticity and fluid dynamics; we choose here the latter approach (see next subsection).

### Appendix A.2.3. Interpretation in Terms of Velocity, Vorticity, and Strain Rates

The heterogeneities of the medium lead to two types of distortions of the flow field, pure strain (without rotation) and pure rotation (without strain). The strain itself can be decomposed into a spherical part (compressive volumetric strain) and an anti-symmetric part (shear strain). In our case, the compressive strain is null. This line of reasoning is borrowed from continuum mechanics and utilizes the fact that the steady-state Darcy flux vector is an incompressible velocity field (divergence-free velocity).

In what follows, we use indicial notations with Einstein's rule of implicit summation on repeated indices. The quantities of interest are the pure strain rate tensor (which is rotation-free):

$$ S_{ij} = \frac{1}{2} \left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right) \qquad (A20) $$

and the pure solid rotation rate tensor (which is strain-free):

$$ R_{ij} = \frac{1}{2} \left( \frac{\partial v_i}{\partial x_j} - \frac{\partial v_j}{\partial x_i} \right) \qquad (A21) $$

Note that the strain and rotation rate tensors are, respectively, symmetric and skew-symmetric. The strain can be further decomposed as the sum of a compressive or spherical strain ($s_{kk}/3\delta_{ij}$), and shear strain or deviatoric strain equal to $s_{ij}$ ($s_{kk}/3\delta_{ij}$). It is easily seen that $(s_{kk}/3) = \text{div}(\mathbf{v})/3 = 0$ due to mass conservation. Thus, the total strain rate equals the shear strain rate for our steady state Darcy flow system without source term (this may not be so in other cases). The total deformation rate $\partial v_i / \partial x_j$ is the sum of (shear) strain rate and (solid) rotation rate:

$$ \frac{\partial v_i}{\partial x_j} = S_{ij} + R_{ij} \qquad (A22) $$

Finally, we relate the vorticity vector given by Equation (A17) to the rotation rate tensor, by the transformations:

$$ \Omega_i = \epsilon_{ijk} R_{kj}; \quad R_{ij} = -\frac{1}{2} \epsilon_{ijk} \Omega_k \qquad (A23) $$

where $\epsilon_{ijk}$ is known as the alternating tensor. It is defined as follows: $\epsilon_{ijk}$ is 0 if any two indices are equal, +1 if $\{i,j,k\}$ is a permutation of $\{1,2,3\}$ in cyclic order, and -1 if $\{i,j,k\}$ is a permutation of $\{1,2,3\}$ in anti-cyclic order. In the three-dimensional case, this yields:

$$ R_{ij} = \frac{1}{2} \begin{bmatrix} 0 & +\Omega_3 & -\Omega_2 \\ -\Omega_3 & 0 & +\Omega_1 \\ +\Omega_2 & -\Omega_1 & 0 \end{bmatrix} \qquad (A24) $$
In the two-dimensional case, the rotation tensor $R_{ij}$ is obtained by elimination of the last row and column of the above matrix. Additionally, observe the factor 1/2 in the last equation (some texts in fluid dynamics define vorticity as 1/2 the curl of velocity).

Now, inserting the rotation rate from Equations (A23) and (A24) in the flux-vorticity Equation (A18) yields the velocity vector equation in the form:

$$ \nabla^2 v_i = 2 \frac{\partial R_{ij}}{\partial x_j} \qquad (A25) $$

Furthermore, inserting the vorticity Equation (A17) into the definition of rotation rate Equations (A23) and (A24) yields a relation between rotation rate, hydraulic conductivity, and velocity (due to Darcy's law) as follows:

$$ R_{ij} = \frac{1}{2} \left( v_i \frac{\partial F}{\partial x_j} - v_j \frac{\partial F}{\partial x_i} \right) \qquad (A26) $$

Finally, inserting this expression in the right-hand side of Equation (A25) gives the desired expanded form of the flux-governing equation:

$$ \nabla^2 v_i = \frac{\partial}{\partial x_j} \left( v_i \frac{\partial F}{\partial x_j} - v_j \frac{\partial F}{\partial x_i} \right) \qquad (A27) $$

where $F$ is the spatially variable log-conductivity $Ln(K)$.

Note that if $K$ is homogeneous (constant in space) then the rotation rate vanishes, the vorticity vanishes, and the flow governing equation becomes $\nabla^2 v_i = 0$.

Additionally, it should be emphasized that the general velocity Equation (A27) is linear, not only with respect to the dependent variable ($v_i$), but also with respect to the heterogeneous log-conductivity field $F = Ln(K)$, whereas the Darcy equation itself is linear in terms of $K$ and not in terms of $F = Ln(K)$.

### Appendix A.3. Analytical Solution with Linearly Varying Head

In this appendix, we derive the analytical solution for the head distribution $H(x_1, x_2)$ within a homogeneous anisotropic domain $\Omega$ under linearly varying heads boundary conditions.

The problem to be solved is:

$$ \nabla(\mathbf{K} \nabla H) = 0 \qquad (A28) $$

or:

$$ K_{11} \frac{\partial^2 H}{\partial x_1^2} + K_{22} \frac{\partial^2 H}{\partial x_2^2} + (K_{12} + K_{21}) \frac{\partial^2 H}{\partial x_2 \partial x_1} = 0 \qquad (A29) $$

with boundary conditions on $\Gamma$:

$$ H(\mathbf{x}) = h_0 - \mathbf{J} \cdot \mathbf{x} \qquad (A30) $$

The solution for the head distribution inside the domain $H(\mathbf{x})$ is simply equal to the linearly varying head prescribed on the boundary whatever is the shape of the domain. It is the solution that both respect the partial differential equation and the boundary conditions. Therefore, in this case the specific discharge is constant in the whole domain:

$$ \mathbf{v}(\mathbf{x}) = \mathbf{K} \mathbf{J} \qquad (A31) $$

### Appendix A.4. Net Surface Flux in 2D

First, we treat the Darcy velocity vector, denoted $\mathbf{A}$ here, as being constant along the domain boundary $\Gamma$. We assume then that only surface fluxes (scalar product $\mathbf{A} \cdot \mathbf{n}$) on the boundary of the domain are available to measurement, and we assume that the
geometry of the domain is known. We want to find a way to estimate the vector **A** from the measurements and the geometry.

Knowing the geometry, we can orient the fluxes depending on the orientation of the boundaries and sum of these fluxes. We define then two components:

$$F_1 = \int_{\Gamma} (\mathbf{A} \cdot \mathbf{n}) (\mathbf{u}_1 \cdot \mathbf{n}) d\gamma \quad (A32)$$

and

$$F_2 = \int_{\Gamma} (\mathbf{A} \cdot \mathbf{n}) (\mathbf{u}_2 \cdot \mathbf{n}) d\gamma \quad (A33)$$

We then introduce the relation

$$\mathbf{A} = A_1 \mathbf{u}_1 + A_2 \mathbf{u}_2 \quad (A34)$$

within the definition of $F_1$ and $F_2$. We illustrate this operation on $F_1$,

$$F_1 = \int_{\Gamma} ((A_1 \mathbf{u}_1 + A_2 \mathbf{u}_2) \cdot \mathbf{n}) (\mathbf{u}_1 \cdot \mathbf{n}) d\gamma \quad (A35)$$

or

$$F_1 = A_1 \int_{\Gamma} (\mathbf{u}_1 \cdot \mathbf{n})^2 d\gamma + A_2 \int_{\Gamma} (\mathbf{u}_2 \cdot \mathbf{n}) (\mathbf{u}_1 \cdot \mathbf{n}) d\gamma \quad (A36)$$

The next step is to define the following integrals, which depend only on the geometry of the domain,

$$N_1 = \int_{\Gamma} (\mathbf{u}_1 \cdot \mathbf{n})^2 d\gamma \quad (A37)$$

$$N_2 = \int_{\Gamma} (\mathbf{u}_2 \cdot \mathbf{n})^2 d\gamma \quad (A38)$$

$$N_{12} = \int_{\Gamma} (\mathbf{u}_1 \cdot \mathbf{n}) (\mathbf{u}_2 \cdot \mathbf{n}) d\gamma \quad (A39)$$

The fluxes $F_1$ and $F_2$ can then be written as:

$$F_1 = A_1 N_1 + A_2 N_{12} \quad (A40)$$

and

$$F_2 = A_1 N_{12} + A_2 N_2 \quad (A41)$$

Solving this linear system of two equations for $A_1$ and $A_2$ gives

$$A_1 = \frac{N_2 F_1 - N_{12} F_2}{N_1 N_2 - N_{12}^2} \quad (A42)$$

and

$$A_2 = \frac{N_1 F_2 - N_{12} F_1}{N_1 N_2 - N_{12}^2} \quad (A43)$$

We see then that if the second members of the two previous equations are used to define the component of the net surface flux vector, then the definition is consistent, i.e., it will systematically provide, for any geometry, the components of a constant vector in the domain from the surface fluxes.

### Appendix A.5. NSF Criteria in 3D

In three dimension, the derivation of the NSF criteria is rather tedious; we must define six geometrical integrals and three oriented net fluxes,

$$N_i = \int_{\Gamma} (\mathbf{u}_i \cdot \mathbf{n})^2 d\gamma; i \in \{1, 2, 3\} \quad (A44)$$
$$N_{ij} = \int_{\Gamma} (\mathbf{u}_{i} \bullet \mathbf{n}) (\mathbf{u}_{j} \bullet \mathbf{n}) d\gamma; i, j \in \{1, 23\} \qquad (A45)$$

$$F_i = \int_{\Gamma} ((\mathbf{v} \bullet \mathbf{n}) (\mathbf{u}_i \bullet \mathbf{n}) d\gamma \qquad (A46)$$

Using all these definitions, the first component of net surface flux is given by:

$$A_1 = \frac{(N_{13}F_2 - N_{12}F_3)(N_{13}N_{12} - N_1N_{23}) - (N_{13}F_1 - N_1F_3)(N_{13}N_2 - N_{12}N_{23})}{(N_{13}N_2 - N_{12}N_{23})(N_{13}^2 - N_1N_3) - (N_{13}N_{12} - N_1N_{23})(N_{13}N_{23} - N_{12}N_3)} \qquad (A47)$$

The other components can be derived following the same principle.

## Appendix A.6. Properties of the Anti-Symmetric Part of a Tensor

Any matrix **K** can be decomposed in the sum of a symmetric **K**<sub>S</sub> and an anti-symmetric **K**<sub>A</sub> matrices

$$\mathbf{K} = \mathbf{K}_S + \mathbf{K}_A \qquad (A48)$$

such that

$$\mathbf{K}_S = \mathbf{K}_S^T \mathbf{K}_S = \mathbf{K}_S^T \text{ and } \mathbf{K}_A = -\mathbf{K}_A^T \qquad (A49)$$

where the superscript T stands for transpose. According to the previous relations, we have directly:

$$\mathbf{K}_S = \frac{\mathbf{K} + \mathbf{K}^T}{2} \qquad (A50)$$

$$\mathbf{K}_A = \frac{\mathbf{K} - \mathbf{K}^T}{2} \qquad (A51)$$

The anti-symmetric matrix is of the form

$$\mathbf{K}_A = \begin{pmatrix} 0 & K_A \\ -K_A & 0 \end{pmatrix} \qquad (A52)$$

A property of an anti-symmetric tensor is that it is invariant by a rotation of the system of coordinate,

$$\mathbf{R}\mathbf{K}_A\mathbf{R}^T = \mathbf{K}_A \qquad (A53)$$

The consequence of this property is that if a tensor of hydraulic conductivity shows an antisymmetric component, it cannot vanish when rotating the system of coordinates. In other words, there will not exist any system of coordinates for which the flow can be collinear with the gradient.

## Appendix A.7. Volume Integral to Surface Integrals

The TDP, VAF, and VAG criteria are defined as volume integrals, while the NSF and VSF criteria are surface integrals. To compare them, it is useful to transform the volume integrals in surface integrals. This is the aim of this section.

### Appendix A.7.1. Green's Gradient and Divergence Formulae

The basic tool is the Green formula, or integration by parts. It can be written in a scalar form as:

$$\int_{\Omega} \frac{\partial f}{\partial x_i} d\omega = \int_{\Gamma} f(\mathbf{x}) n_i(\mathbf{x}) d\gamma \qquad (A54)$$

where $n_i(\mathbf{x})$, the ith component of the normal vector to the elementary surface $d\gamma$. The same expression can be written in vectorial form as:

$$\int_{\Omega} \nabla f d\omega = \int_{\Gamma} f(\mathbf{x}) \mathbf{n} d\gamma \qquad (A55)$$

with **n**, the normal vector to the elementary surface $d\gamma$.
The above identity (“gradient integral”) can be derived as follows from the classical flux-divergence theorem (which is another Green identity). Let us start therefore with the flux divergence theorem:

$$ \int_{\Omega} \nabla \mathbf{W} d\omega = \int_{\Gamma} \mathbf{W} \mathbf{n} d\gamma \qquad (A56) $$

Usually, **W** represents a flux or velocity vector, and the left-hand term represents the volume integral of the “flux divergence” *div*(*W*), while the right-hand side term is the net flux through the boundary of the domain. Now, the trick is to replace the vector **W** by a diagonal matrix **W**: (we use here exceptionally double underlining for matrices):

$$ \underline{\underline{\mathbf{W}}} = \begin{bmatrix} f & & \\ & f & \\ & & f \end{bmatrix} \qquad (A57) $$

Observe now that we have the following identities (using here vector and indicial notations):

$$ \mathrm{div}(\underline{\underline{\mathbf{W}}}) = \nabla \mathbf{W} = \nabla^T \mathbf{W} = \nabla^T \begin{bmatrix} f & & \\ & f & \\ & & f \end{bmatrix} \Rightarrow \qquad (A58) $$

$$ \mathrm{div}(\underline{\underline{\mathbf{W}}}) = \begin{bmatrix} \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \end{bmatrix} \begin{bmatrix} f & & \\ & f & \\ & & f \end{bmatrix} = \begin{bmatrix} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \\ \frac{\partial f}{\partial z} \end{bmatrix} $$

This leads finally to the desired “gradient” theorem, or generalized integration by parts, given above in Equations (A54) and (A55).

## Appendix A.7.2. Simplification of the VAG

In the general case, we obtain a simplification (or another formulation) of the VAG through a direct application of the Green formula (Sánchez-Vila et al., 1995):

$$ \langle \nabla h \rangle_{\mathrm{VAG}} = \frac{1}{\Omega} \int_{\Omega} \nabla h \, d\omega = \frac{1}{\Omega} \int_{\Gamma} h(\mathbf{x}) \, n \, d\gamma \qquad (A59) $$

### Specialization for Linear Head Conditions (“Immersion”)

However, we note also, more specifically, that an additional simplification occurs if the boundary conditions are of linear head distributions (“immersion” conditions). We can write them as

$$ h(\mathbf{x}) = -\mathbf{J} \cdot \mathbf{x} + h_0 \mathbf{x} \in \Gamma \qquad (A60) $$

Using Equations (A59) and (A60), the volume average head gradient in the domain becomes:

$$ \langle \nabla h \rangle_{\mathrm{VAG}} = \frac{1}{\Omega} \int_{\Gamma} (-\mathbf{J} \cdot \mathbf{x} + h_0) \mathbf{n} \, d\gamma \qquad (A61) $$

or

$$ \langle \nabla h \rangle_{\mathrm{VAG}} = \frac{1}{\Omega} \int_{\Gamma} (-\mathbf{J} \cdot \mathbf{x}) \mathbf{n} \, d\gamma + \frac{1}{\Omega} \int_{\Gamma} h_0 \mathbf{n} \, d\gamma \qquad (A62) $$

The second term vanishes because $h_0$ is constant. We can then apply Green’s formula back to transform the first surface integral into a volume integral, and we obtain:

$$ \langle \nabla h \rangle_{\mathrm{VAG}} = \frac{1}{\Omega} \int_{\Omega} -\mathbf{J} \, d\omega = -\mathbf{J} \qquad (A63) $$

To sum up: the volume average hydraulic gradient inside the domain, $-\langle \nabla h \rangle_{\mathrm{VAG}}$, is equal to the constant head gradient **J** used to define the prescribed head boundary conditions on the boundary surface of the block. The result is independent of the geometry
of the domain and of the spatial heterogeneities within the block. The proof may appear
strange; we started with a volume integral to transform it in a surface integral, and finally
we transformed it back into a volume integral. The key is that the internal distribution
of the heads $h(x)$ is unknown for the microscale heterogeneous medium. Only the head
distribution on the boundary is known (prescribed). The first transform allows to replace
the unknown heads inside the domain by the heads on the boundary. When these boundary
heads are known (prescribed), we can move back to the interior of the domain with a known
function which is different from the actual heads in the domain.

Specialization for Rectangular Geometry

When the geometry of the domain is a parallelepiped, the surface integral can be
separated on every face of the domain. To illustrate this principle, the integral can be split
into four terms in the case of the rectangular body previously introduced. First, let us
consider only one component of the average gradient:

$$
\langle \nabla h \rangle_1 = \frac{1}{\Omega} \int_{\Gamma} h \mathbf{u}_1 \cdot \mathbf{n} \, d\gamma \qquad (A64)
$$

The two integrals on the faces $\Gamma_{A2}$ and $\Gamma_{B2}$ vanish, as the normal vectors are perpendicular to the unit vector $\mathbf{u}_1$. Finally, the $x_1$ component of the gradient is in this case:

$$
\langle \nabla h \rangle_1 = \frac{\frac{1}{\Gamma_{B1}} \int_{\Gamma_{B1}} h \, d\gamma - \frac{1}{\Gamma_{A1}} \int_{\Gamma_{A1}} h \, d\gamma}{L_1} \quad (A65)
$$

Specialization for Periodic Conditions

When the boundary conditions are of periodic type, then the average head gradient within the domain is equal to the imposed head gradient at the boundary of the domain. In the rectangular example, we have:

$$
\langle \nabla h \rangle_{VAG} = \frac{1}{\Omega} \int_{\Gamma} h(\mathbf{x}) \, \mathbf{n} \, d\gamma = - \left( \frac{\frac{\Delta h_1}{L_1}}{\frac{\Delta h_2}{L_2}} \right) = - \sum_{i=1,2} \left( \frac{\Delta h_i}{L_i} \right) \mathbf{u}_i \quad (A66)
$$

Specialization for Permeametric Conditions

When the boundary conditions are of permeameter type, the average head gradient in
the direction of imposed flow is simply the one prescribed by the boundary conditions be-
tween the two opposite faces divided by the separating distance between the two opposite
faces. Assuming, for example, a gradient in the $x_1$ direction, we have:

$$
\langle \nabla h \rangle_1 = - \frac{\Delta h_1}{L_1} \qquad (A67)
$$

In the perpendicular direction(s), the average head gradient is not imposed by the
boundary conditions, but it is also different from zero in general. For example, in the $x_2$
direction, we have:

$$
\langle \nabla h \rangle_2 = \frac{\frac{1}{\Gamma_{B2}} \int_{\Gamma_{B2}} h \, d\gamma - \frac{1}{\Gamma_{A2}} \int_{\Gamma_{A2}} h \, d\gamma}{L_2} \quad (A68)
$$

Polyhedra with N Faces at Constant Head

In the case of a polyhedral domain that is limited by *n* planar boundary faces, we have

$$
\langle \nabla h \rangle_{VAG} = \frac{\sum_{i=1}^{n} h_i \mathbf{n}_i S_i}{\Omega} \qquad (A69)
$$

with $h_i$ the prescribed head on face $i$, $\mathbf{n}_i$ the normal vector and $S_i$ the area of face $i$.
## Appendix A.7.3. Simplification of the VAF Expression

The volume average flux is simplified in two steps. First, we consider the $x_i$ component of the average specific discharge flux (Darcy velocity):

$$ \langle v \rangle_i = \left( \frac{1}{\Omega} \int_{\Omega} \mathbf{v} \, d\omega \right) \bullet \mathbf{u}_i = \frac{1}{\Omega} \int_{\Omega} \mathbf{v} \bullet \mathbf{u}_i \, d\omega \quad (A70) $$

Then, an integration by part (Green's formula) leads to:

$$ \langle v \rangle_i = \frac{1}{\Omega} \left( \int_{\Gamma} x_i \mathbf{v} \bullet \mathbf{n} d\gamma - \int_{\Omega} \mathbf{u}_i \nabla \cdot \mathbf{v} d\omega \right) \quad (A71) $$

Since the flow is steady-state $\nabla \cdot \mathbf{v} = 0$, we get

$$ \langle v \rangle_i = \frac{1}{\Omega} \int_{\Gamma} x_i (\mathbf{v} \bullet \mathbf{n}) \, d\gamma \quad (A72) $$

or

$$ \langle \mathbf{v} \rangle_{\text{VAF}} = \frac{1}{\Omega} \int_{\Gamma} \mathbf{x} (\mathbf{v} \bullet \mathbf{n}) \, d\gamma \quad (A73) $$

## Specialization for Rectangular Geometry

If we assume the rectangular shape, and if we note $q_i$ the flux across the face $i$ and $x_{Gi}$ the center of gravity of the face $i$ weighted by the fluxes on the face, we have:

$$ q_i = \int_{\Gamma_i} \mathbf{v} \bullet \mathbf{n} \, d\gamma \quad (A74) $$

and

$$ x_{Gi} = \frac{\int_{\Gamma_i} \mathbf{x} \mathbf{v} \bullet \mathbf{n} \, d\gamma}{q_i} \quad (A75) $$

We can then express the average specific discharge as follows (as in Sánchez-Vila et al., 1995):

$$ \langle v \rangle_1 = \frac{1}{\Omega} \int_{\Gamma} x_i \mathbf{v} \bullet \mathbf{n} \, d\gamma = \frac{q_{A1}x_1^{A1} + q_{B1}x_1^{B1} + q_{A2}x_{G1}^{A2} + q_{B2}x_{G1}^{B2}}{\Omega} \quad (A76) $$

With $x_1^{A1}$ being the (constant) $x_1$ coordinate of face A1, and $x_{G1}^{A2}$ the $x_1$ coordinate of the center of gravity of face A2 weighted by the fluxes (see Figure 1 in the text). Note that in the case of linearly varying head boundary conditions, these equations cannot be simplified further. See below for the case of permeametric conditions.

## Specialization for Permeametric Conditions

When the boundary conditions are of permeameter type, then the previous expression can be further simplified because

$$ q_{A2} = q_{B2} = 0, \quad q_{A1} = -q_{B1}, \quad x_{1G}^{A2} = x_{1G}^{B2} = 0 \quad (A77) $$

Using these relations, we obtain the following expressions:

$$ \langle v \rangle_1 = \frac{q_{B1}}{\Gamma_1} \quad (A78) $$

and

$$ \langle v \rangle_2 = q_{B1} \frac{(x_{G2}^{B1} - x_{G2}^{A1})}{\Gamma_1} \quad (A79) $$
Specialization for Periodic Conditions

When the boundary conditions are periodic, the distribution of the fluxes along opposite surfaces are identical. Consequently, we have the following identities:

$$q_{A1} = -q_{B1}, \quad x_{G2}^{A1} = x_{G2}^{B1}, \quad q_{A2} = -q_{B2}, \quad x_{G1}^{A2} = x_{G1}^{B2} \qquad (A80)$$

and therefore, the volume average fluxes are simply:

$$\langle v \rangle_1 = \frac{q_{B1}}{\Gamma_1} \qquad (A81)$$

and

$$\langle v \rangle_2 = \frac{q_{B2}}{\Gamma_2} \qquad (A82)$$

Appendix A.7.4. Simplification of the Total Dissipated Power (TDP)

The total dissipated power is also simplified thanks to an integration by parts (Boe, 1994, Njifenjou, 1994). Let us first recall the definition of P:

$$P = -\rho g \int_{\Omega} (\mathbf{v} \cdot \nabla h) d\omega \qquad (A83)$$

Integration by part (or Green's identity) gives:

$$P = -\rho g \left[ \int_{\Gamma} h (\mathbf{v} \cdot \mathbf{n}) d\gamma - \int_{\Omega} h (\nabla \cdot \mathbf{v}) d\omega \right] \qquad (A84)$$

Since the velocity field is divergence-free ($\nabla \mathbf{v} = 0$) we obtain:

$$P = -\rho g \int_{\Gamma} h \mathbf{v} \cdot \mathbf{n} d\gamma \qquad (A85)$$

Rectangular Geometry

Let us assume again a rectangular geometry. The integral can be decomposed on all the elemental surfaces:

$$-\frac{P}{\rho g} = \int_{\Gamma_{A1}} h(\mathbf{v} \cdot \mathbf{n}) d\gamma + \int_{\Gamma_{B1}} h(\mathbf{v} \cdot \mathbf{n}) d\gamma + \int_{\Gamma_{A2}} h(\mathbf{v} \cdot \mathbf{n}) d\gamma + \int_{\Gamma_{B2}} h(\mathbf{v} \cdot \mathbf{n}) d\gamma \qquad (A86)$$

Specialization for Permeametric Conditions

When the boundary conditions are of permeameter type, the head is constant on two opposite faces and the normal fluxes vanishes on the other boundary faces ("lateral" faces). We obtain, then:

$$-\frac{P}{\rho g} = h_0 q_{A1} + (h_0 - \Delta h_1) q_{B1} \qquad (A87)$$

Since $q_{A1} = -q_{B1}$, this yields:

$$P = \rho g \Delta h_1 q_{B1} \qquad (A88)$$

Specialization for Periodic Conditions

When the boundary conditions are periodic, the simplification requires to regroup the integrals into two subgroups or "components", related to the opposite faces. For example, if we take the pair of opposite faces ($A_1, B_1$), we define the "component" $P_1$ of P as follows
(Note: the name “component” is only used for convenience here; P is a scalar; and the «components» (P₁, P₂) are such that P = P₁ + P₂.):

$$ -\frac{P_1}{\rho g} = \int_{\Gamma_{A1}} h(\mathbf{v}\cdot\mathbf{n}) d\gamma + \int_{\Gamma_{B1}} h(\mathbf{v}\cdot\mathbf{n}) d\gamma \quad (A89) $$

These integrals are expressed as a function of the coordinate system.

$$ -\frac{P_1}{\rho g} = \int_{x_2=0}^{x_2=L_1} h(0, x_2) \mathbf{v}(0, x_2) \cdot (-\mathbf{u}_1) dx_2 + \int_{x_2=0}^{x_2=L_1} h(L_1, x_2) \mathbf{v}(L_1, x_2) \cdot \mathbf{u}_1 dx_2 \quad (A90) $$

We can group the two integrals and because of the boundary conditions v(0, x₂) = v(L₁, x₂), we obtain:

$$ -\frac{P_1}{\rho g} = \int_{x_2=0}^{x_2=L_1} [h(L_1, x_2) - h(0, x_2)] \mathbf{v}(0, x_2) \cdot \mathbf{u}_1 dx_2 \quad (A91) $$

or

$$ P_1 = \rho g \Delta h_1 q_{B1} \quad (A92) $$

The same applies for the second direction, and therefore, the total dissipated power is obtained from P = P₁ + P₂, which yields finally:

$$ P = \rho g \Delta h_1 q_{B1} + \rho g \Delta h_2 q_{B2} \quad (A93) $$

Specialization for Linearly Varying Head Boundary Conditions (Gradient Conditions)

When the boundary conditions are linearly varying heads, the hydraulic head on Γ is

$$ h(\mathbf{x}) = h_0 - \mathbf{J} \cdot \mathbf{x} \quad (A94) $$

The total dissipated power is then

$$ P = -\rho g \int_{\Gamma} (h_0 - \mathbf{J} \cdot \mathbf{x}) (\mathbf{v} \cdot \mathbf{n}) d\gamma \quad (A95) $$

or

$$ P = \rho g \mathbf{J} \cdot \int_{\Gamma} \mathbf{x} (\mathbf{v} \cdot \mathbf{n}) d\gamma - \rho g h_0 \int_{\Gamma} (\mathbf{v} \cdot \mathbf{n}) d\gamma \quad (A96) $$

The last integral vanishes because ∇v = 0 within Ω. We recognize that the first integral is the volume average flux within the domain (see Equation (A73)). Therefore, we obtain:

$$ P = -\rho g \mathbf{J} \cdot \langle \mathbf{v} \rangle_{VAF} \quad (A97) $$

## Appendix A.8. Dual Permeability/Resistivity Upscaling Procedure

This Appendix presents more details on the dual permeability/resistivity homogenization procedure, as a complement to the main text (Section 4.3).

The dual permeability/resistivity was developed in particular for Darcy flow in randomly heterogeneous spatially correlated media on infinite domains, based on spectral perturbation solutions of the permeability and resistivity-based Darcy flow equations, along with an exponential extrapolation of the resulting macro-permeability and macro-resistivity [75,103]. In particular, two criteria are required: positivity of the resulting effective permeability and resistivity tensors (K and R must be positive definite); and reciprocity relation (R must be the inverse of K).

It should be noted that this dual upscaling procedure is applicable to both single-phase and two-phase flows in locally isotropic porous media (with conductivity and resistivity dependence on capillary pressure in the case of two-phase flow); in addition, it is potentially applicable for evaluating the equivalent permeability of finite size blocks without necessarily assuming random-type heterogeneity. The key is matching the (approximate) macro-permeability with the (approximate) macro-resistivity.
Our presentation below illustrates the method in particular for statistical continua in infinite domains.

First, the local microscale flow equations are cast in both permeability and resistivity form:

$$ \mathbf{v} = -k(\mathbf{x}) \nabla h \operatorname{div}(\mathbf{v}) = 0 \qquad (A98a) $$

$$ \nabla h = -r(\mathbf{x}) \mathbf{v} \operatorname{div}(\mathbf{v}) = 0 \qquad (A98b) $$

where the scalar resistivity $r(\mathbf{x})$ is defined as the inverse of the scalar conductivity:

$$ r(\mathbf{x}) = 1/k(\mathbf{x}) $$

Note that these Darcy relations, Equation (A98a) for permeability, as well as Equation (A98b) for resistivity, involve the two vectors $\mathbf{v}$ and $\nabla h$. Another approach is possible in terms solely of the velocity vector $\mathbf{v}$ by using the identity $\nabla \times \nabla h = 0$ (see Section 2.2.2 and Appendix A.2), but this “velocity approach” is not directly used here for permeability/resistivity matching.

Next, an upscaling procedure is applied to both equations in Equation (A98a,b). In the case at hand, the log-coefficients $F = Ln k(\mathbf{x}) = -Ln r(\mathbf{x})$ are used, and the upscaling procedure is applied twice, $\mathbf{v} = -k(\mathbf{x}) \nabla h$ and to $\nabla h = -r(\mathbf{x}) \mathbf{v}$. Upscaling is based here on expressing ensemble averages of $\mathbf{v}(\mathbf{x})$ and $\nabla h(\mathbf{x})$ based on perturbation approximations using the standard deviation $\sigma = \sigma_{Ln k}$ as a “small parameter”. The quantity $Ln k(\mathbf{x}) = -Ln r(\mathbf{x})$ is a statistically homogeneous random field, with well-defined auto-correlation function and spectral density.

If the coordinate system coincides with the principal axes of statistical anisotropy of the random medium, the resulting effective permeability and effective resistivity obtained from this procedure, to second order accuracy in terms of $\sigma$, are two distinct tensors of the form:

$$ k_{ij}^{\mathrm{EFF}(2)} \approx k_G \left\{ \delta_{ij} + \frac{\sigma^2}{2} X_{ij}^* \right\} \rightarrow r_{ij}^{\mathrm{EFF}(2)} \approx \frac{1}{k_G} \left\{ \delta_{ij} - \frac{\sigma^2}{2} X_{ij}^* \right\} \qquad (A99) $$

where the label (2) in $k_{ij}^{\mathrm{EFF}(2)}$ and $r_{ij}^{\mathrm{EFF}(2)}$ indicate second order perturbation approximations, the star superscript “*” indicates principal values, and $X_{ij}$ is a microgeometric tensor that does not depend on the hydraulic gradient (for single phase flow), but only on statistical geometry and correlation structure of the random porous medium (in the spectral perturbation theory it is a Fourier integral of the spectral density function of the random field $F = Ln k(\mathbf{x}) = -Ln r(\mathbf{x})$). Note that $X_{ij}^*$ (with the star*) represents the components of tensor $X_{ij}$ expressed in the coordinate system aligned with the principal axes of the porous medium (this is indicated by the * superscript). Thus, in Equation (A99), we have taken the coordinate system to be aligned with the principal axes of statistical anisotropy. In that case, the effective $k_{ij}^{\mathrm{EFF}(2)}$ and $r_{ij}^{\mathrm{EFF}(2)}$ are diagonal, and it is easily seen that (a) the principal values $r_{ij}^{\mathrm{EFF}(2)}$ are not the inverse of $k_{ij}^{\mathrm{EFF}(2)}$, and (b) one of them at least could become negative for large enough $\sigma$. Thus, at this point, remarkably, the reciprocity relation and the positivity criterion are not satisfied.

Next, to improve robustness towards larger variability (larger $\sigma$), the above approximations are modified according to the following “exponential extrapolation” (a conjectural extrapolation of perturbation results, presented in Ababou [102], Gelhar [32], and proved to be exact in the 1D case in Ababou [111]):

$$ k_{ij}^{\mathrm{EFF(EX)}} \approx k_G \exp\left\{\frac{\sigma^2}{2} X_{ij}^*\right\} \rightarrow r_{ij}^{\mathrm{EFF(EX)}} \approx \frac{1}{k_G} \exp\left\{-\frac{\sigma^2}{2} X_{ij}^*\right\} \qquad (A100) $$

However, the above relation holds only for $X_{ij}^*$ expressed in the principal coordinate system. Due to the exponential extrapolation, the tensorial nature of the effective coefficients
is now lost: the rule to transform the permeability and resistivity to a new coordinate system not aligned with principal axes ... is not known. This is because the exponential $\exp\{X_{ij}\}$ of a tensor like $X_{ij}$ is not a true tensor even if $X_{ij}$ is a tensor.

Therefore, a last modification is suggested in order to conserve the tensorial nature of effective coefficients using the following empirical formulation:

$$k_{ij}^{\text{EFF(EX.TENSOR)}} \approx k_G \mathbf{P}^{\text{T}} \exp\left\{\frac{\sigma^2}{2} X_{ij}^*\right\} \mathbf{P} \quad (\text{A101})$$

$$r_{ij}^{\text{EFF(EX.TENSOR)}} \approx \frac{1}{k_G} \mathbf{P}^{\text{T}} \exp\left\{-\frac{\sigma^2}{2} X_{ij}^*\right\} \mathbf{P}$$

This last formulation preserves the tensorial nature of effective permeability and resistivity, that is, it satisfies the tensorial transformation rules of the form $T \rightarrow T' = P^T T P$. The permeability and resistivity tensors are now expressed in any coordinate system, as a transformation from the (starred*) principal system, using the passage matrix $\mathbf{P}$ (see Section 2.4.1 and Appendix A.1.1).

To sum up, this last formulation (Equation (A101)) achieves several desirable properties: (i) it is thought to ensure the robustness of the small parameter approximation thanks to the exponential extrapolation; (ii) it satisfies both the reciprocity relation and the positivity criterion for permeability and resistivity; and (iii) it preserves the tensorial nature of effective permeability and resistivity.

## Nomenclature

<table><thead><tr><th>Abbreviation</th><th>Description</th></tr></thead><tbody><tr><td>NSF</td><td>Net Surface Flux</td></tr><tr><td>VSF</td><td>Vectorial Surface Flux</td></tr><tr><td>VAF</td><td>Volume Average Flux</td></tr><tr><td>VAG</td><td>Volume Average Gradient</td></tr><tr><td>TDP</td><td>Total Dissipated Power</td></tr><tr><td>PRM</td><td>PeRMeameter boundary conditions</td></tr><tr><td>GRD</td><td>Constant head GRaDient immersion boundary conditions</td></tr><tr><td>FLX</td><td>Constant FLuX immersion boundary conditions</td></tr><tr><td>PRD</td><td>PeRioDic boundary conditions</td></tr><tr><td>SKN</td><td>Skin immersion</td></tr><tr><td>GRD'</td><td>Piecewise constant GRaDient conditions</td></tr><tr><td>FLX'</td><td>Piecewise constant FLuX conditions</td></tr><tr><td>PRM'</td><td>The "dual permeameter" (non-zero fluxes on two opposite faces, linear heads on the other two faces).</td></tr><tr><th>Symbol</th><th>Description</th></tr><tr><td>h, H</td><td>Local and global hydraulic head (m)</td></tr><tr><td>j, J</td><td>Local and global hydraulic gradient -∇h, -∇H</td></tr><tr><td>k, K</td><td>Local and global hydraulic conductivity tensor (m/s)</td></tr><tr><td>P</td><td>Total dissipated power (see note below) (J/s).</td></tr><tr><td>v, V</td><td>Local and global Darcy velocity, or specific discharge rate (m/s)</td></tr></tbody></table>

## References

1. Baas, J.H.; Hailwood, E.A.; McCaffrey, W.D.; Kay, M.; Jones, R. Directional petrological characterisation of deep-marine sandstones using grain fabric and permeability anisotropy: Methodologies, theory, application and suggestions for integration. *Earth-Sci. Rev.* **2007**, *82*, 101–142. [CrossRef]

2. Burger, R.L.; Belitz, K. Measurement of anisotropic hydraulic conductivity in unconsolidated sands: A case study from a shoreface deposit, Oyster, Virginia. *Water Resour. Res.* **1997**, *33*, 1515–1522. [CrossRef]

3. Dewhurst, D.N.; Brown, K.M.; Clennell, M.B.; Westbrook, G.K. A comparison of the fabric and permeability anisotropy of consolidated and sheared silty clay. *Eng. Geol.* **1996**, *42*, 253–267. [CrossRef]

4. Greenkorn, R.A.; Johnson, C.; Shallenberger, L.K. Directional permeability of heterogeneous anisotropic porous media. *Soc. Pet. Eng. J.* **1964**, *4*, 124–132. [CrossRef]
5. Grimestad, G. Analysis of data from pumping tests in anisotropic aquifers—Equations and graphical solutions. *Water Resour. Res.* **1995**, *31*, 933–941. [CrossRef]

6. Hurst, A.; Rosvoll, K.J. Permeability variations in sandstones and their relationship to sedimentary structures. In *Reservoir Characterization II*; Lake, L.W., Carroll, H.B., Wesson, T.C., Eds.; Academic Press: Cambridge, MA, USA, 1991; pp. 166–196.

7. Hutta, J.J.; Griffiths, J.C. Directional permeability of sandstones; a test of technique. *Prod. Mon.* **1955**, *19*, 12–24.

8. Moore, P.J. Determination of permeability anisotropy in a two-way permeameter. *Geotech. Test. J.* **1979**, *2*, 167–169.

9. Motz, L.H. Multiple-pumped-well aquifer test to determine the anisotropic properties of a karst limestone aquifer in Pasco County, Florida, USA. *Hydrogeol. J.* **2009**, *17*, 855–869. [CrossRef]

10. Papadopulos, I.S. Nonsteady flow to a well in an infinite anisotropic aquifer. *Hydrol. Fract. Rocks* **1965**, *21*–31.

11. Phonate, V.; Malik, R.S.; Kumar, S. Anisotropy in some alluvial soils. *Ann. Arid. Zone* **2001**, *40*, 425–427.

12. Rayne, T.W.; Mickelson, D.M. Sediment fabric and anisotropy of hydraulic conductivity in Sandy till, Wisconsin, USA. *Hydrol. Res.* **1996**, *27*, 161–174. [CrossRef]

13. van den Berg, E.H.; de Vries, J.J. Influence of grain fabric and lamination on the anisotropy of hydraulic conductivity in unconsolidated dune sands. *J. Hydrol.* **2003**, *283*, 244–266. [CrossRef]

14. Borghi, A.; Renard, P.; Courrioux, G. Generation of 3D spatially variable anisotropy for groundwater flow simulations. *Groundwater* **2015**, *53*, 955–958. [CrossRef] [PubMed]

15. Medici, G.; Smeraglia, L.; Torabi, A.; Botter, C. Review of modeling approaches to groundwater flow in deformed carbonate aquifers. *Groundwater* **2021**, *59*, 334–351. [CrossRef]

16. Bennett, J.P.; Haslauer, C.P.; Cirpka, O.A. The impact of sedimentary anisotropy on solute mixing in stacked scour-pool structures. *Water Resour. Res.* **2017**, *53*, 2813–2832. [CrossRef]

17. Sánchez-Vila, X.; Carrera, J. Directional effects on convergent flow tracer tests. *Math. Geol.* **1997**, *29*, 551–569. [CrossRef]

18. Michael, H.A.; Voss, C.I. Controls on groundwater flow in the Bengal Basin of India and Bangladesh: Regional modeling analysis. *Hydrogeology* **2009**, *17*, 1561–1577. [CrossRef]

19. Hemker, K.; van den Berg, E.; Bakker, M. Ground water whirls. *Ground Water* **2004**, *42*, 234–242. [CrossRef]

20. Stauffer, F. Impact of highly permeable sediment units with inclined bedding on solute transport in aquifers. *Adv. Water Res.* **2007**, *30*, 2194–2201. [CrossRef]

21. Zhao, P.; Xie, L.; He, B.; Liu, J. Anisotropic permeability influencing the performance of free CH4 and free CO2 during the process of CO2 sequestration and enhanced gas recovery (CS-EGR) from shale. *ACS Sustain. Chem. Eng.* **2021**, *9*, 914–926. [CrossRef]

22. Weijermars, R.; Khanal, A. High-resolution streamline models of flow in fractured porous media using discrete fractures: Implications for upscaling of permeability anisotropy. *Earth-Sci. Rev.* **2019**, *194*, 399–448. [CrossRef]

23. Bernabé, Y. On the measurement of permeability in anisotropic rocks. In *Fault Mechanics and Transport Properties of Rocks*; Evans, B., Wong, T.-F., Eds.; Academic Press: Cambridge, MA, USA, 1992; pp. 147–167.

24. Bieber, M.T.; Rasolofosaon, P.; Zinszner, B.; Zamora, M. Measurement and overall characterization of permeability anisotropy by tracer injection. *Rev. L'institut Français Pétrole* **1996**, *51*, 333–347. [CrossRef]

25. Renard, P.; Genty, A.; Stauffer, F. Laboratory determination of the full tensor of permeability of a sample. *J. Geophys. Res.* **2001**, *106*, 26443–26452. [CrossRef]

26. Clavaud, J.-B.; Maineult, A.; Zamora, M.; Rasolofosaon, P.; Schlitter, C. Permeability anisotropy and its relations with porous medium structure. *J. Geophys. Res. Earth Surf.* **2008**, *113*, B01202. [CrossRef]

27. Hsieh, P.; Neuman, S.; Stiles, G.K.; Simpson, E.S. Field determination of the three-dimensional hydraulic conductivity tensor of anisotropic media 1. Methodology and application to fractured rocks. *Water Resour. Res.* **1985**, *21*, 1667–1676. [CrossRef]

28. Moench, A.F. Flow to a well of finite diameter in a homogeneous, anisotropic water table aquifer. *Water Resour. Res.* **1997**, *33*, 1397–1407. [CrossRef]

29. Cihan, A.; Zhou, Q.; Birkholzer, J.T.; Kraemer, S.R. Flow in horizontally anisotropic multilayered aquifer systems with leaky wells and aquitards. *Water Resour. Res.* **2014**, *50*, 741–747. [CrossRef]

30. Matheron, G. Éléments Pour Une Théorie des Milieux Poreux; Masson & Cie: Paris, France, 1967.

31. Dagan, G. *Flow and Transport in Porous Formations*; Springer Science & Business Media: Berlin, Germany, 1989.

32. Gelhar, L.W. *Stochastic Subsurface Hydrology*; Prentice Hall: Englewood Cliffs, NJ, USA, 1993; p. 390.

33. Wen, X.-H.; Gómez-Hernández, J.J. Upscaling hydraulic conductivities in heterogeneous media: An overview. *J. Hydrol.* **1996**, *183*, ix–xxxii. [CrossRef]

34. Renard, P.; de Marsily, G. Calculating equivalent permeability: A review. *Adv. Water Res.* **1997**, *20*, 253–278. [CrossRef]

35. Farmer, C.L. Upscaling: A review. *Int. J. Numer. Methods Fluids* **2002**, *40*, 63–78. [CrossRef]

36. Sanchez-Vila, X.; Guadagnini, A.; Carrera, J. Representative hydraulic conductivities in saturated groundwater flow. *Rev. Geophys.* **2006**, *44*, RG3002. [CrossRef]

37. Durlofsky, L.J. Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media. *Water Resour. Res.* **1991**, *27*, 699–708. [CrossRef]

38. Kfoury, M.; Ababou, R.; Noetinger, B.; Quintard, M. Upscaling Fractured Heterogeneous Media: Permeability and Mass Exchange Coefficient. *J. Appl. Mech.* **2006**, *73*, 41–46. [CrossRef]

39. Lang, P.S.; Paluszny, A.; Zimmerman, R.W. Permeability tensor of three-dimensional fractured porous rock and a comparison to trace map predictions. *J. Geophys. Res. Solid Earth* **2014**, *119*, 6288–6307. [CrossRef]
40. Long, J.C.; Remer, J.S.; Wilson, C.R.; Witherspoon, P.A. Porous media equivalents for networks of discontinuous fractures. *Water Resour. Res.* **1982**, *18*, 645–658. [CrossRef]

41. Massier, S. *Simulation Directe D'Ecoulements Sur des Réseaux Bidimensionnels de Fractures: Conductivité Equivalente*; Stage de 2ème année, Ecole Centrale de Lyon: Lyon, France, 1994.

42. Pickup, G.E.; Ringrose, P.S.; Jensen, J.L.; Sorbie, K.S. Permeability Tensors for Sedimentary Structures. *Math. Geol.* **1994**, *26*, 227–250. [CrossRef]

43. Pouya, A.; Fouche, O. Permeability of 3D discontinuity networks: New tensors from boundary-conditioned homogenisation. *Adv. Water Res.* **2009**, *32*, 303–314. [CrossRef]

44. Scandelli, H.; Ahmadi-Senichault, A.; Levet, C.; Lachaud, J. Computation of the Permeability Tensor of Non-Periodic Anisotropic Porous Media from 3D Images. *Transp. Porous Media* **2022**, *142*, 669–697. [CrossRef]

45. Afra, S.; Gildin, E. Permeability Parametrization Using High Order Singular Value Decomposition (HOSVD). In Proceedings of the 12th International Conference on Machine Learning and Applications, Miami, FL, USA, 4–7 December 2013; pp. 188–193. [CrossRef]

46. Moslehi, M.; de Barros, F.P.J.; Ebrahimi, F.; Sahimi, M. Upscaling of solute transport in disordered porous media by wavelet transformations. *Adv. Water Resour.* **2016**, *96*, 180–189. [CrossRef]

47. Amaziane, B.; Koebbe, J. JHomogenizer: A computational tool for upscaling permeability for flow in heterogeneous porous media. *Comput. Geosci.* **2006**, *10*, 343–359. [CrossRef]

48. Dagan, G.; Fiori, A.; Jankovic, I. Upscaling of flow in heterogeneous porous formations: Critical examination and issues of principle. *Adv. Water Resour.* **2013**, *51*, 67–85. [CrossRef]

49. Flodin, E.A.; Durlofsky, L.J.; Aydin, A. Upscaled models of flow and transport in faulted sandstone: Boundary condition effects and explicit fracture modelling. *Pet. Geosci.* **2004**, *10*, 173–181. [CrossRef]

50. Giudici, M.; Vassena, C. About the symmetry of the upscaled equivalent transmissivity tensor. *Math. Geol.* **2007**, *39*, 399–408. [CrossRef]

51. Pouya, A.; Courtois, A. Definition of the permeability of fractured rock masses by homogenisation methods. *Comptes Rendus Geosci.* **2002**, *334*, 975–979. [CrossRef]

52. Zhang, Y.; Liu, B.; Gable, C.W. Homogenization of Hydraulic Conductivity for Hierarchical Sedimentary Deposits at Multiple Scales. *Transp. Porous Med.* **2011**, *87*, 717–737. [CrossRef]

53. Neuman, S.P.; Orr, S. Prediction of steady state flow in nonuniform geologic media by conditional moments: Exact nonlocal formalism, effective conductivities and weak approximation. *Water Resour. Res.* **1993**, *29*, 341–364. [CrossRef]

54. Stam, J.M.T.; Zijl, W. Modeling permeability in imperfectly layered media. II. A two-dimensional application of block-scale permeability. *Math. Geol.* **1992**, *24*, 885–904.

55. Wen, X.-H.; Durlofsky, L.J.; Edwards, M.G. Use of border regions for improved permeability upscaling. *Math. Geol.* **2003**, *35*, 521–547. [CrossRef]

56. Darcel, C.; Bour, O.; Davy, P.; De Dreuzy, J.R. Connectivity properties of two-dimensional fracture networks with stochastic fractal correlation. *Water Resour. Res.* **2003**, *39*, 1272. [CrossRef]

57. De Dreuzy, J.-R.; De Boiry, P.; Pichot, G.; Davy, P. Use of power averaging for quantifying the influence of structure organization on permeability upscaling in on-lattice networks under mean parallel flow. *Water Resour. Res.* **2010**, *46*, W08519. [CrossRef]

58. Hestir, K.; Long, J.C.S. Analytical expressions for the permeability of random 2-dimensional Poisson fracture networks based on regular lattice percolation and equivalent media theories. *J. Geophys. Res. Solid* **1990**, *95*, 21565–21581. [CrossRef]

59. Kiraly, L. Anisotropie et hétérogénéité de la perméabilité dans les calcaires fissurés. *Eclogae Geol. Helv.* **1969**, *62*, 613–619.

60. Snow, D.T. Anisotropic permeability of fractured media. *Water Resour. Res.* **1969**, *5*, 1273–1289. [CrossRef]

61. Ababou, R.; Renard, P. Equivalent permeability tensor in fractured media: An algebraic approach. In Proceedings of the MAMERN11: 4th International Conference on Approximation Methods and Numerical Modelling in Environment and Natural Resources, Saidia, Morocco, 23–26 May 2011; Amaziane, B., Barrera, D., Mraoui, H., Rodriguez, M.L., Sbibih, D., Eds.; Universidad de Granada: Granada, Spain, 2011; p. 4, ISBN 078-84-338-5230-4.

62. Renard, P.; Ababou, R. Relation between the definition and properties of the equivalent permeability tensor in heterogeneous and fractured porous media. In Proceedings of the MAMERN09 3rd International Conference on Approximation Methods & Numerical Modelling in Environment & Natural Resources, Pau, France, 8–11 June 2009; Amaziane, B., Barrera, D., Sbibih, D., Eds.; Editorial Universidad de Granada: Granada, Spain, 2009; p. 4, ISBN 978-84-338-5006-5.

63. Cañamón, I.; Rajeh, T.; Ababou, R.; Marcoux, M. Topological analysis of 3D fracture networks: Graph representation and percolation threshold. *Comput. Geotech.* **2022**, *142*, 104556. [CrossRef]

64. Rajeh, T.; Ababou, R.; Marcoux, M.; Cañamon, I. Fast upscaling of the hydraulic conductivity of three-dimensional fractured porous rock for reservoir modeling. *Math. Geosci.* **2019**, *51*, 1037–1074. [CrossRef]

65. Barker, A.T.; Lee, C.S.; Vassilevski, P.S. Spectral upscaling for graph Laplacian problems with application to reservoir simulation. *SIAM J. Sci. Comput.* **2017**, *39*, S323–S346. [CrossRef]

66. Ferrandon, J. Les lois de l'écoulement de filtration. *Génie Civ.* **1948**, *125*, 24–28.

67. Scheidegger, A.E. Directional permeability of porous media to homogeneous fluids. *Geofis. Pura E Appl.* **1954**, *28*, 75–90. [CrossRef]

68. Ababou, R. Three-Dimensional Flow in Random Porous Media. Ph.D. Thesis, Ralph Parsons Laboratory, Massachusetts Institute of Technology, Cambridge, MA, USA, 1988. Volume 2. p. 833.
69. Zijl, W.; Nawalany, M. Natural Groundwater Flow. Lewis Publishers: Boca Raton, FL, USA, 1993; p. 321.

70. Akpoji, A.G.; Ababou, R.; De Smedt, F. Stochastic flow analysis of flux variances in heterogeneous porous media. In *Computational Methods in Water Resources X, Tenth International Conference held at Universität Heidelberg, Germany, July 1994*; Peters, A., Wittum, G., Herrling, B., Meissner, U., Brebbia, C.A., Gray, W.G., Pinder, G.F., Eds.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1994; Volume 1, pp. 391–398. ISBN 0-7923-2935-X.

71. Diestel, R. *Graph Theory*, 2nd ed.; Springer: Cham, Switzerland, 2000.

72. Wilson, R.J. *Introduction to Graph Theory*, 4th ed.; Pearson Education: Essex, UK, 1996.

73. Shvidler, M.I. *Filtration Flows in Heterogeneous Media*; A Statistical Approach; Consultants Bureau: New York, NY, USA, 1964.

74. Quintard, M.; Whitaker, S. Ecoulement monophasique en milieu poreux: Effets des hétérogénéités locales. *J. Mécanique Théorique Appliquée* **1987**, *6*, 691–726.

75. Fadili, A.; Ababou, R. Dual homogenization of immiscible steady two-phase flows in random porous media. *Water Resour. Res.* **2004**, *40*. [CrossRef]

76. Aris, R. *Vector, Tensors, and the Basic Equations of Fluid Mechanics*; Dover: Mineola, NY, USA, 1989; 1st edition by Prentice Hall: Englewood Cliffs, NJ, USA, 1962.

77. Bear, J. *Dynamics of Fluids in Porous Media*; Dover Publication Inc.: Mineola, NY, USA, 1988; [Republished from American Elsevier, New York, 1972].

78. Bamberger, A. *Approximation des Coefficients D'Opérateurs Elliptiques, Stable Pour la G-Convergence*; Ecole Polytechnique: Paris, France, 1977.

79. Gallouët, T.; Guérillot, D. Averaged heterogeneous porous media by minimisation of the error on the flow rate. In Proceedings of the ECMOR IV-4th European Conference on the Mathematics of Oil Recovery, Røros, Norway, 7–10 June 1994.

80. Renard, P. Modélisation des Ecoulements en Milieux Poreux Hétérogènes, Calcul des Perméabilités Equivalentes. Ph.D. Thesis, Ecole des Mines de Paris, Paris, France, 1997. Mémoires des Sciences de la Terre, Volume 32.

81. Rajeh, T.; Ababou, R.; Marcoux, M. Conductivité hydraulique équivalente d'une roche poreuse fracturée dans un réservoir géothermique profond. In Proceedings of the 13èmes Journées d'études des Milieux Poreux, Anglet, France, 12–14 October 2016. (hal-01394520).

82. Cañamón, I. *Analysis and Modeling of Coupled Thermo-Hydro-Mechanical Phenomena in 3D Fractured Media*; Institut National Polytechnique: Toulouse, France, 2006.

83. Cañamón, I.; Ababou, R.; Elorza, F.J. A 3-dimensional homogenized model of coupled thermo-hydro-mechanics for nuclear waste disposal in geologic media. In Proceedings of the European Nuclear Conference 2007 (ENC) on Nuclear Waste Modeling, Brussels, Belgium, 16–20 September 2007.

84. Ababou, R.; Valera, I.C.; Poutrel, A. Macro-permeability distribution and anisotropy in a 3D fissured and fractured clay rock: 'Excavation Damaged Zone' around a cylindrical drift in Callovo-Oxfordian Argilite (Bure). *Phys. Chem. Earth, Parts A/B/C* **2011**, *36*, 1932–1948. [CrossRef]

85. Pouya, A. Equivalent permeability tensors of a finite heterogeneous block. *Comptes Rendus Geosci.* **2005**, *337*, 581–588. [CrossRef]

86. Ababou, R.; Bagtzoglou, A.C. BIGFLOW: A Numerical Code for Simulating Flow in Variably Saturated, Heterogeneous Geologic Media (*Theory and User's Manual, Version 1.1*), Report NUREG/CR-6028; U.S. Nuclear Regulatory Commission, Government Printing Office: Washington, DC, USA, 1993; p. 139. Available online: http://www.osti.gov/bridge/servlets/purl/10168217-yoTsuT/10168217.pdf (accessed on 11 May 2022).

87. Gómez-Hernández, J.J.; Journel, A.G. Stochastic characterization of grid-block permeabilities: From point values to block tensors. *SPE Form. Eval.* **1994**, *9*, 83–99. [CrossRef]

88. Philip, J.R. Issues in flow and transport in heterogeneous porous media. *Transp. Porous Media* **1986**, *1*, 319–338. [CrossRef]

89. Kapoor, V.; Kitanidis, P. Dilution of nonreactive solutes in heterogeneous porous media. In *Scale Invariance and Scale Dependence in Hydrology*; Sposito, G., Ed.; Cambridge University Press: New York, NY, USA, 1996.

90. Cardwell, W.T.; Parsons, R.L. Average permeabilities of heterogeneous oil sands. *Trans. AIME* **1942**, *160*, 34–42. [CrossRef]

91. Warren, J.E.; Price, H.S. Flow in heterogeneous porous media. *Soc. Pet. Eng. J.* **1961**, *1*, 153–169. [CrossRef]

92. Rubin, Y.; Gómez-Hernández, J.J. A stochastic approach to the problem of upscaling of conductivity in disordered media: Theory and unconditional numerical simulations. *Water Resour. Res.* **1990**, *22*, 691–701. [CrossRef]

93. Sánchez-Vila, X.; Girardi, J.P.; Carrera, J. A synthesis of approaches to upscaling of hydraulic conductivities. *Water Resour. Res.* **1995**, *31*, 867–882. [CrossRef]

94. Loc'h, G.L. An efficient strategy for combining the permeabilities: Practical applications on a simulated reservoir. In *Geostatistics*; Armstrong, M., Ed.; Springer: Dordrecht, The Netherlands, 1989; Volume 2, pp. 557–568.

95. Indelman, P. Upscaling of permeability of anisotropic heterogeneous formations, 3, applications. *Water Resour. Res.* **1993**, *29*, 935–943. [CrossRef]

96. Bøe, Ø. Analysis of an upscaling method based on conservation of dissipation. *Transp. Porous Media* **1994**, *17*, 77–86. [CrossRef]

97. Ababou, R.; Renard, P. Exact algebraic approach to flow calculation and permeability upscaling in 2D fracture networks using graph operators. Oral Communication in Session 51 on Fractured geological media and fracture networks: Flow, graphs, morphology. In Proceedings of the IAMG 2022, 21st Annual Conference of the International Association for Mathematical Geosciences, Nancy, France, 29 August–3 September 2022.

98. Strang, G. A framework for equilibrium equations. *SIAM Rev.* **1988**, *30*, 283–297. [CrossRef]
99. Strang, G. *Computational Science and Engineering*; Wellesley-Cambridge Press: Wellesley, MA, USA, 2007.

100. Indelman, P.; Dagan, G. Upscaling of heterogeneous formations: General approach and application to isotropic media. *Transp. Porous Media* **1993**, *12*, 61–183. [CrossRef]

101. Njifenjou, A. Expression en termes d'énergie pour la perméabilité absolue effective. *Rev. De L'institut Français Du Pétrole* **1994**, *49*, 345–358. [CrossRef]

102. Ababou, R. Stochastic Homogenization and Effective Conductivity Tensor in Random Anisotropic Media. In Proceedings of the ICIAM'91 2nd International Conference on Industrial & Applied Mathematics—Minisymposium M146 on Darcy Flow in Composite and Random Porous Media SIAM, Washington, DC, USA, 8–12 July 1991; pp. 300–302.

103. Ababou, R. Stochastic velocity-pressure fields in random porous media: Statistical expansions, solutions, and homogenization. In Proceedings of the International Conference on Industrial and Applied Mathematics ICIAM'99, Edinburgh, UK, 5–9 July 1999.

104. Cañamón, I.; Elorza, F.J.; Ababou, R. 3D Fracture Networks: Optimal Identification and Reconstruction. In Proceedings of the IAMG'06: International Association of Mathematical Geology, XIth International Congress, Université de Liège, Liège, Belgium, 3–8 September 2006.

105. Kasap, E.; Lake, L.W. Calculating the effective permeability tensor of a gridblock. *SPE Form. Eval.* **1990**, *5*, 192–200. [CrossRef]

106. Ababou, R. Random porous media flow on large 3-D grids: Numerics, performance and application to homogenization. In *IMA Volumes in Mathematics and its Applications Environmental Studies: Mathematical, Computational and Statistical Analysis*; Wheeler, M.F., Ed.; Springer Science & Business Media: Berlin, Germany, 1996; pp. 1–25.

107. Householder, A.S. *The Theory of Matrices in Numerical Analysis*. Dover Books on Mathematics; Dover Publication Inc.: Mineola, NY, USA, 2006.

108. Christie, M.A.; Blunt, M.J. Tenth SPE comparative solution project: A comparison of upscaling techniques. *SPE Reserv. Eval. Eng.* **2001**, *4*, 308–317. [CrossRef]

109. Noetinger, B.; Jarrige, N. A quasi steady state method for solving transient Darcy flow in complex 3D fractured networks. *J. Comput. Phys.* **2012**, *231*, 23. [CrossRef]

110. Gradshtein, I.S.; Rhyzik, I.M. *Tables of Integrals, Series, and Products*, 7th ed.; Jeffrey, A., Zwillinger, D., Eds.; Academic Press: Cambridge, MA, USA, 2007; p. 1171.

111. Ababou, R. Solution of Stochastic Groundwater Flow by Infinite Series, and Convergence of the One-Dimensional Expansion. *Stoch. Hydrol. Hydraul.* **1994**, *8*, 139–155. [CrossRef]