# Equivalent Permeability Tensor of Heterogeneous Media: Upscaling Methods and Criteria (Review and Analyses)

Philippe Renard, Rachid Ababou

► To cite this version:

Philippe Renard, Rachid Ababou. Equivalent Permeability Tensor of Heterogeneous Media: Upscaling Methods and Criteria (Review and Analyses). Geosciences, 2022, 12 (7), pp.269. 10.3390/geosciences12070269. [https://doi.org/10.3390/geosciences12070269](https://doi.org/10.3390/geosciences12070269). [hal-04483557](https://wwwHAL.fr/hal-04483557)

HAL Id: hal-04483557

https://ut3-toulouseinp.hal.science/hal-04483557v1

Submitted on 29 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Review

# Equivalent Permeability Tensor of Heterogeneous Media: Upscaling Methods and Criteria (Review and Analyses)

Philippe Renard<sup>1,2,*</sup> and Rachid Ababou<sup>3</sup>

**Abstract:** When conducting numerical upscaling, either for a fractured or a porous medium, it is important to account for anisotropy because in general, the resulting upscaled conductivity is anisotropic. Measurements made at different scales also demonstrate the existence of anisotropy of hydraulic conductivity. At the "microscopic" scale, the anisotropy results from the preferential flatness of grains, presence of shale, or variation of grain size in successive laminations. At a larger scale, the anisotropy results from preferential orientation of highly conductive geological features (channels, fracture families) or alternations of high and low conductive features (stratification, bed-ding, crossbedding). Previous surveys of homogenization techniques demonstrate that a wide variety of approaches exists to define and calculate the equivalent conductivity tensor. Consequently, the resulting equivalent conductivities obtained by these different methods are not necessarily equal, and they do not have the same mathematical properties (some are symmetric, others are not, for example). We present an overview of different techniques allowing a quantitative evaluation of the anisotropic equivalent conductivity for heterogeneous porous media, via numerical simulations and, in some cases, analytical approaches. New approaches to equivalent permeability are proposed for heterogeneous media, as well as discontinuous (composite) media, and also some extensions to 2D fractured networks. One of the main focuses of the paper is to explore the relations between these various definitions and the resulting properties of the anisotropic equivalent conductivity, such as tensorial or non-tensorial behavior of the anisotropic conductivity; symmetry and positiveness of the conductivity tensor (or not); dual conductivity/resistivity tensors; continuity and robustness of equivalent conductivity with respect to domain geometry and boundary conditions. In this paper, we emphasize some of the implications of the different approaches for the resulting equivalent permeabilities.

**Keywords:** upscaling; permeability tensor; porous media; fractured media; anisotropy; Darcy's law; heterogeneity

## 1. Introduction

The hydraulic conductivity of geological media is often anisotropic as it has been shown in laboratory and field experiments on a wide variety of scales and a wide variety of porous or fractured rocks [1–13]. At the microscopic (pore) scale, the anisotropy is due to the preferential flatness of grains, the presence of shale, or the variation of grain size in successive laminations. At a larger scale, the anisotropy results from preferential orientation of geological features such as fracture families, stratifications, or channels [14,15]. Anisotropy is a key parameter controlling the behavior of solutes such as contaminants in porous materials [16,17]. An example demonstrating the crucial importance of anisotropy for regional hydrogeology studies is the study of Arsenic contamination in the Bengal basin [18]. Anisotropy can trigger the existence of groundwater whirls [19,20] that may

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Renard, P.; Ababou, R. Equivalent Permeability Tensor of Heterogeneous Media: Upscaling Methods and Criteria (Review and Analyses). *Geosciences* **2022**, *12*, 269. https://doi.org/10.3390/geosciences12070269

Academic Editors: Helder I. Chaminé and Jesus Martinez-Frias

Received: 12 May 2022

Accepted: 20 June 2022

Published: 1 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Centre for Hydrogeology and Geothermics, University of Neuchâtel, Rue Emile-Argand 11, 2000 Neuchâtel, Switzerland</footnote>

<footnote n="2">² Department of Geosciences, Sem Sælands vei, University of Oslo, 1 Geologibygningen, 0371 Oslo, Norway</footnote>

<footnote n="3">³ Institut de Mécanique des Fluides de Toulouse, CNRS—University of Toulouse, 2 Allée Camille Soula, 31400 Toulouse, France; ababou@imft.fr</footnote>

<footnote>* Correspondence: philippe.renard@unine.ch</footnote>
have a large influence on reactive transport and dilution processes. Anisotropy can im-
pact the efficiency of greenhouse gases (such as CO2), sequestration in deep geological
formations [21], and heat extraction by geothermal boreholes [22].

Therefore, there is a need for methods and tools allowing to estimate anisotropy in
the laboratory [23–26] or in the field [27–29], or to evaluate anisotropy through upscaling
techniques (e.g., when a detailed geological model has been built but is too detailed for
direct numerical modeling). When anisotropy is put in this perspective, one realizes
that some common techniques can be used to interpret permeametric experiments in the
laboratory as well as numerical experiments based on geological models.

Those techniques have been essentially developed over the last 50 years in the frame-
work of upscaling, either on finite domains or finite “blocks” (equivalent block conductivity,
as in the present work) or in infinite domain (effective hydraulic conductivity). Analyses
and reviews of these methods can be found, for example, in [30–36]. Their applications
to estimate tensorial permeabilities in numerical experiments or in the laboratory have
been studied by many authors (including, for example, [23,25,37–44]). Other works have
focused on optimal estimation of heterogeneous reservoir properties such as permeability
k(x, y). For instance, techniques such as truncated Proper Orthogonal Decomposition
(POD) have been used in this context, as well as High Order Singular Value Decomposi-
tions (HOSVD), which lead typically to a coarsened optimal estimation of permeability
distribution (e.g., [45]). Other orthogonal decompositions such as wavelet decompositions
have been used for intermediate upscaling of permeability k(x, y) with application to solute
transport upscaling (e.g., [46]).

Note that within the upscaling literature, one can distinguish *effective properties* and *equivalent block properties*. The effective properties emerge when the size of the heterogeneities is much smaller than the size of the sample and when there is some statistical homogeneity within the sample or some geometrical periodicity. When these conditions are not met, the equivalent conductivity is not any more an intrinsic property of the medium, and to avoid confusion with effective permeability one uses, in general, terminology such as *equivalent block permeability*, or simply *block permeability*. In the following, we will consider only block permeability tensors. The two concepts converge when the conditions of emergence of an effective permeability are met.

The block hydraulic conductivity, or block permeability tensor, obtained from an experiment (either numerical or physical) can be very different depending on the technique used to analyze the experiment (type of averaging, for example) and on the type of experiment itself (essentially the type of boundary conditions imposed on the sample). This is especially the case if the porous sample exhibits anisotropic heterogeneity features having a scale comparable to the size of the sample. This effect is described in most of the reviews cited above, and in several other papers [39,43,47–52].

In particular, it is known that depending on the technique used, one can obtain non-
symmetrical tensors (e.g., [53–55]). This has raised controversial issues. For example, if
the hydraulic conductivity tensor is not symmetric, it implies that there is no direction
such that the Darcy velocity is parallel to the hydraulic gradient, and one cannot identify
anymore the usual principal directions of anisotropy. Note, however, that it is possible to
analyze a non-symmetric permeability tensor in terms of two sets of principal directions
and principal components.

Similarly, it is possible to obtain tensors that may not be positive definite. For instance,
the permeability tensor may be indefinite; or else it may be semi-definite positive (i.e., pos-
itive in some directions but zero in some other directions). In the most counter-intuitive
situations, one can obtain a permeability tensor having one negative eigenvalue, that is, a
negative principal permeability (see example in Section 6.2). This implies that the mean
flow along this principal direction will go “upgradient”, from zones of low hydraulic head
to zones of high hydraulic head. In addition, if the permeability tensor is not definite posi-
tive, it implies that there may be some directions of flow such that the energy dissipated by
the viscous forces in the medium can become zero or even negative (a situation that has no
physical sense). Indeed, the only way to ensure that the dissipated power is always strictly positive when the head gradient is not equal to zero is to ensure that the permeability tensor is definite positive.

To better understand the relations between the definitions and properties of the permeability tensors obtained from upscaling procedures or averaging techniques in numerical or laboratory experiments, we present an overview of the most frequent techniques, and we discuss new approaches as well. Specifically, the main focus of this work is to explore the relations between these various definitions and the resulting properties of the anisotropic equivalent conductivity.

Fracture networks constitute a special case of interest [40,51,56–60]. If the porous matrix is assumed impervious, then flow takes place only through an intersecting network of discrete conductive objects (the percolating subset of the fracture network). If the network is considered to be 2D, then flow occurs through a bond network of links and nodes (or its percolating subset). An algebraic approach to the equivalent permeability of such a 2D network was proposed by Ababou and Renard [61,62]. The case of 3D planar fractures is more complex. The percolating subset and other properties of a 3D network of planar disc fractures were analyzed recently by Cañamón et al. [63]. In these works, the rock matrix was assumed impervious, and the emphasis was on the topological graph properties of the fracture networks. Additionally, because of this extreme type of discontinuity (flowing fractures vs. impervious matrix), it is often necessary to use special dedicated formulations for flux averages and equivalent conductivities. The intermediate case of a 3D permeable rock matrix traversed by planar disc fractures was also treated by Rajeh, Ababou et al. [64]; their study of permeability upscaling took into account both the permeability contrast Kf/Km and the density of the 3D fracture set. They also considered the asymptotic case Kf/Km → ∞, and they finally expressed the upscaled tensorial permeability of the matrix/fracture medium in terms of a critical density of fractures (the so-called critical “exclusion density” based on the exclusion volume of fractures). Other works that focus on fractured media include Long [40], Pouya and Fouché [43], and Lang et al. [39]. Barker et al. [65] consider the spectral graph properties of networks and developed a reservoir simulation application involving a Finite Volume discretization of a 3D heterogeneous continuum, re-interpreted as a 3D network of links.

In the following, we consider three types of heterogeneous media:

*   *continuous porous media:* this includes the special cases of uniform (homogenous) media, non-uniform (inhomogeneous) media with continuously varying properties, and spatially correlated random media;

*   *composite porous media:* this includes various multilayered media, as well as more general Boolean porous media (binary random mixtures, thresholded random fields, etc.) where permeability is non-zero everywhere;

*   *fractured networks:* we include in this category the special case of 2D fractured media with impervious porous matrix; this case is equivalent to the problem of flow on a network.

In the remainder of this work, we focus to a larger extent on the first two types of heterogeneous media mentioned above: (i) continuous porous media, and (ii) composite porous media, and to a lesser extent on the third case (iii) of fractured media/2D networks.

The remainder of this paper is organized as follows:

In the next section (Section 2), we present our assumptions of flow and geometry, and we show how the conductivity tensor can be computed from the results of numerical experiments. This requires that one defines both the boundary conditions (experimental setup) and the equivalence criteria to be used to define the equivalent conductivity, usually through some kind of averaging. In Section 3, we present some special analytical solutions for Darcy flow in finite domains, which provide insights on the effects of permeability variability on the flow pattern and can sometimes lead to analyses of the corresponding equivalent permeability tensor. In Section 4, various average quantities and equivalence criteria are proposed to define and calculate the equivalent hydraulic conductivity of
heterogeneous blocks. Combining these average quantities and criteria with various types of boundary conditions or flow experiments leads to several different ways of calculating the equivalent permeability tensor. In Section 6, numerical experiments are implemented on 10 samples of heterogeneous porous media, including continuously variable random field conductivity patterns, as well as composite or random Boolean conductivity patterns. The special case of 2D fracture networks is treated separately, analytically, in Section 6.2. Section 7 presents analytical and algebraic proofs concerning the properties of equivalent conductivity tensors obtained by some of the previous methods, and the relations between the different methods. A brief conclusive section is provided in Section 8. It is followed by several Appendices, which complete some of the developments presented in the text concerning heterogeneous flow patterns, averages, boundary conditions, etc.

**2. Definitions and Assumptions on Flow and Geometry**

Defining upscaling techniques and characterizing their properties requires first defining the following items: the local flow model; the global or upscaled flow model; a formal or general definition of the upscaled permeability tensor; a series of (possibly alternative) equivalence criteria or homogenization criteria; and a description of the boundary conditions used to calculate the equivalent permeability or conductivity tensor. Note: a brief introduction to tensorial permeability will be presented in Section 2.4, with more details in Appendix A.1, concerning second rank tensors in general, then conductivity tensors in particular, ending up with the concept of directional conductivity (when flow is governed by the tensorial form of Darcy's law).

**2.1. Flow Domain: Block Geometry**

We do not fix a priori the geometry of the domain of interest. It can either be a simple geometrical shape (rectangle, parallelepiped, triangle, tetrahedron, etc.) or completely general. In some situations, we may have to use a parallelepipedic geometry. In such a case, we will use the conventions introduced in Figure 1.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Geometrical sketch of: (a) a general 2D domain $\Omega$, with *n* the outgoing normal vector with respect to the boundary $\Gamma$. (b) Rectangular domain geometry and size. (c) Faces and normal vectors in a rectangular geometry. The letter *A* stands for the faces closest to the origin. The letter *B* corresponds to the opposite face. The index (1,2, or 3) stands for the axis perpendicular to the face. (d) The same principle is used to denote the faces of a 3D parallelepipedic domain as well.</figcaption></figure>
## 2.2. Local Flow Models

At the macroscopic scale, the distribution of permeability within the porous medium (or the spatial distribution of the fractures and their apertures) is assumed to be known. The geometry of the domain is also assumed to be known. Note: in the case of a fractured porous medium, we would consider two local flow models, one for the porous medium and one for the fractures; however, here we only provide indications for the case of 2D fracture networks neglecting the permeability of the porous matrix.

### 2.2.1. Head-Based Flow Equations for a Porous Medium

The porous medium is assumed to be water-saturated (or liquid saturated more generally). The flow is assumed to follow Darcy's law:

$$ \mathbf{v}(\mathbf{x}) = -\mathbf{k}(\mathbf{x}) \nabla h(\mathbf{x}) \qquad (1) $$

where $\mathbf{v}(\mathbf{x})$ [LT⁻¹] is the local Darcy velocity vector or specific discharge (often denoted $\mathbf{q}(\mathbf{x})$); $h$ [L] is the local total hydraulic head, and $\mathbf{k}(\mathbf{x})$ [LT⁻¹] is the local hydraulic conductivity tensor, which is assumed symmetric and positive-definite [66,67].

$$ \mathbf{k} = \begin{pmatrix} k^{11} & k^{12} & k^{13} \\ k^{12} & k^{22} & k^{23} \\ k^{13} & k^{23} & k^{33} \end{pmatrix} \qquad (2) $$

Note that the local, “hydraulic gradient” is often defined as $\mathbf{j}(\mathbf{x}) = -\nabla h(\mathbf{x})$ in the literature, so that Darcy's law can also be expressed as $\mathbf{v}(\mathbf{x}) = \mathbf{k}(\mathbf{x}) \mathbf{j}(\mathbf{x})$.

In addition, the Darcy velocity $\mathbf{v}(\mathbf{x})$ obeys a steady state mass conservation equation (assuming here incompressible fluid):

$$ \mathrm{div}(\mathbf{v}(\mathbf{x})) = \nabla \cdot \mathbf{v}(\mathbf{x}) = 0 \qquad (3) $$

Finally, the local head-based flow equation is obtained by inserting Darcy's law Equation (1) into Equation (3):

$$ - \nabla \cdot (\mathbf{k}(\mathbf{x}) \nabla h(\mathbf{x})) = 0 \qquad (4) $$

In what follows, we will usually assume that the local conductivity is isotropic; therefore $\mathbf{k}(\mathbf{x})$ is a spatially distributed scalar rather than a tensor (locally). (In the remainder of this paper, we use interchangeably the terms “permeability” and “hydraulic conductivity” for convenience, although strictly speaking, k & K represent here local & block scale hydraulic conductivity [m/s].) The local total hydraulic head $h$ [L] is defined as the sum of pressure potential plus gravitational potential, converted into water column height as follows:

$$ h = \frac{p - p_{\mathrm{ATM}}}{\rho g} + z \qquad (5) $$

where $z$ is the vertical elevation (pointing upwards), $p$ is pore water pressure, $p_{\mathrm{ATM}}$ is air pressure (atmospheric), $g$ is the acceleration of gravity, and $\rho$ is the density of liquid water (or the density of the incompressible liquid that saturates the porous medium).

### 2.2.2. Velocity-Based Flow Equations for a Porous Medium

Alternatively, a “flux”-based flow equation governing directly the Darcy velocity in a heterogeneous porous medium has been proposed ([68] Section 4.3 therein, [69] Section 3.5 therein). This equation is obtained by re-stating Darcy's law Equation (1) as follows and assuming from now on that the local permeability is isotropic (scalar field $k$):

$$ \mathbf{v}(\mathbf{x}) / k(\mathbf{x}) = -\nabla h(\mathbf{x}) \qquad (6) $$
We can then infer that the curl of $v(x)/k(x)$ must be null because the curl of a gradient is always null. By definition, the curl (rotational) of a vector field $w(x)$ is the vector product (or cross product) of the vectorial $\nabla$ operator with vector $w$. The velocity equation (which is a direct consequence of Darcy's law) is given by:

$$ \nabla \times \{v(x)/k(x)\} = 0 \qquad (7) $$

Finally, this rotational Equation (7), together with the divergence Equation (3), form a velocity-based system of vector equations, which can be used to solve directly for the velocity field $v(x)$ instead of the head-based Equation (5). Once solved, the hydraulic gradient $-\nabla h(x)$ can then be obtained from Equation (6), and the scalar head field can then be obtained by integration from any boundary point where the head is prescribed.

A more detailed derivation of these velocity or flux-based equations is presented later below in Section 3.1 and in Appendix A.2. In the case of randomly heterogeneous media, the consequences of these local velocity-based equations, as opposed to head-based equations, were studied by Akpoji et al. [70].

The main consequence of this alternative formulation, in the context of equivalent permeability, is that the averaged version of this flow equation may not yield the same equivalent permeability as the head-based flow equations (see discussion in Section 2.3).

### 2.2.3. Head-Based Flow Equations for a 2D Fracture Network (Indications)

In the fractures, the flow is assumed to follow a linear law such as the cubic law. In two dimensions, the fracture network is represented by a graph [71,72] of edges and nodes (as shown in Figure 2).

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Example of a simple 2D fracture network, represented by a directed graph ("digraph"). Black circles such as ● represent the vertices of the graph (i.e., the nodes of the intersecting fracture network), and the line segments with numbers [2] represent the "arcs" or "edges" of the graph (i.e., the fracture segments or links between nodes). The links [ℓ] are interconnected to each other via vertices (nodes). The vertices can be numbered in any order, and similarly for the edge numbers. For a fracture network, the edge directions can also be set up using any predefined convention (see text).</figcaption></figure>

The flux $q_i$ [LT⁻³] along the edge "i" is also denoted $v_i$ for convenience (being similar to a velocity). It is modeled according to either Darcy's law or Poiseuille's law, by:

$$ q_i \text{ or } v_i = c_i(h_j - h_k) \qquad (8) $$

with $c_i$ [L²T⁻¹] the conductance of edge "i", $h_j$ and $h_k$ [L] the hydraulic heads at the nodes corresponding to the extremities of edge "i", $l_i$ the distance between the two end-nodes of edge "i". Note that we have used here the single index "i" to label the edges or links, and the double index $(j, k)$ to label their two end-nods. Alternatively, it is possible to use
the double index $(j, k)$ to label conductivities (e.g., $k_{jk}$) and conductances (e.g., $c_{jk}$) of the fracture kinks $(j, k)$.

The conservation equation is then written for each node of the graph. This law (known as Kirchhoff's law) states that the balance of the fluxes of all the edges $(j, k)$ connected to a given node $(j)$ is equal to zero (zero net flux at any given node):

$$ \text{For each node } (j): \quad \sum_{k \in \mathcal{B}(j)} q_k = 0 \quad \text{or} \quad \sum_{k \in \mathcal{B}(j)} v_k = 0 \qquad (9) $$

where $k \in \mathcal{B}(j)$ represents the indices $k$ of all nodes connected to the node $j$ by a fracture link $(j, k)$. For example, in a regular hexagonal network, each node has six neighbors, and the Kirchhoff law above expresses that, at any node of the network, the sum of all six velocities (ongoing and outgoing) must be zero.

## 2.3. Upscaled Flow Model with Equivalent Conductivity

It is assumed that both in the case of a fractured or a porous medium, the large-scale flow can be described by a macroscopic model analogous to the local Darcy and mass conservation equations, but with macroscale or block scale quantities (we use capital letters for these macroscale quantities):

$$ \mathbf{V}(\mathbf{x}) = -\mathbf{K}(\mathbf{x}) \nabla H(\mathbf{x}) \qquad (10) $$

$$ \mathbf{K} = \begin{pmatrix} K^{11} & K^{12} & K^{13} \\ K^{21} & K^{22} & K^{23} \\ K^{31} & K^{32} & K^{33} \end{pmatrix} \qquad (11) $$

$$ \nabla \mathbf{V}(\mathbf{x}) = 0 \qquad (12) $$

hence, the head-based flow equation is obtained by inserting Darcy's law Equation (10) into Equation (12):

$$ - \nabla \cdot (\mathbf{K}(\mathbf{x}) \nabla H(\mathbf{x})) = 0 \qquad (13) $$

In these macroscale equations, $\mathbf{K}(\mathbf{x})$ [LT$^{-1}$] is the upscaled or equivalent hydraulic conductivity (possibly non-symmetric), $H(\mathbf{x})$ [L] is the macroscale hydraulic head within the upscaled block, and $\mathbf{V}(\mathbf{x})$ is the macroscale specific discharge (Darcy velocity). The macroscale "hydraulic gradient" is defined as $\mathbf{J}(\mathbf{x}) = -\nabla H(\mathbf{x})$.

In all cases, it is assumed that the macroscale $\mathbf{V}$, $\nabla H$, or $\mathbf{J}$ are continuous and smoothly varying functions within the domain. We will assume in many cases that the macroscale quantities $\mathbf{K}$, $\mathbf{V}$, and $\nabla H$ can be considered spatially constant within the block (although this is not necessary), and in such cases, the macroscale head $H(\mathbf{x})$ is a linear function of space. In particular, note that the macroscale velocity $\mathbf{V}$ is either constant or smoothly varying, while in contrast, the local Darcy velocity vector $\mathbf{v}(\mathbf{x})$ can be highly variable and even discontinuous, for instance in composite or fractured media where the tangential component of $\mathbf{v}(\mathbf{x})$ can be discontinuous across a fracture or a material discontinuity within the block. Finally, note also that $\mathbf{K}$, $\mathbf{V}$, $\nabla H$ are not necessarily defined as simple arithmetic averages of the local variables $\mathbf{k}(\mathbf{x})$, $\mathbf{v}(\mathbf{x})$, $\nabla H(\mathbf{x})$.

Several authors, including ([30,73,74], [68] Chapter 4 therein), and others, have demonstrated that such a model governs the large-scale flow if the following conditions are met:

a. the heterogeneities have a small size compared to domain size (or, for a random medium, the spatial correlation scales of the medium are small with respect to domain size, which is analogous to the theoretical "infinite domain" hypothesis); and

b. the flow field is "homogeneous" at the macro-scale or domain scale (or, in the case of a random medium, the flow field should be statistically homogeneous up to second order).

However, when dealing with practical situations, the previous conditions may not be fulfilled. For instance, the mean large-scale flow could be radial, as would occur around a pumping well, or the heterogeneous structures of the medium could be as large as the
size of the sample. If this would be the case, implementing the simplified macroscopic phenomenological law of Equations (10) and (12) with constant **K**, **V**, $\nabla H$ would be mistaken, because we know it is not theoretically exact when the cited conditions on inhomogeneity and/or scales are not met.

As mentioned in Section 2.2.2, the averaged velocity-based flow equations may not yield the same equivalent permeability as the head-based flow equations. 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 is then to match the equivalent permeability (obtained from head-based equations) to equivalent resistivity (obtained from velocity-based equations). We refer to Fadili and Ababou [75] for a similar dual permeability-resistivity upscaling approach in the context of randomly heterogeneous media for single-phase and two-phase flow. A more general matching procedure could be attempted for arbitrarily heterogeneous porous and/or fractured blocks, but this possible extension will not be developed further here.

## 2.4. Tensorial Conductivity & Directional Conductivity Ellipse

This sub-section presents brief definitions of second rank tensors (such as the conductivity tensor $K_{ij}$) and several related concepts: principal axes, principal values, directional "flux" and "gradient" conductivities, anisotropic conductivity ellipses. This topic is relevant for both microscale (local) and macroscale (upscaled) flow models. For references, see [76] concerning vectors and tensors in fluid mechanics, and [77] concerning principal and directional conductivities. A more detailed presentation of the algebraic aspects of this subsection can be found in Appendix A.1.

### 2.4.1. Vectors, Second Rank Tensors, and Tensorial Conductivity Vectors, Tensors, and Transformation Rules

A vector is a first rank tensor. A second rank tensor $\mathbf{T}$ or $T_{ij}$ ($i = 1, \dots, N$; $j = 1, \dots, N$) can be viewed as a linear application and can be represented by a matrix, for instance (in 2D with $N = 2$):

$$ \mathbf{T} \text{ as a square matrix: } \mathbf{T} = \begin{bmatrix} T_{11} & T_{12} \\ T_{21} & T_{22} \end{bmatrix} \qquad (14) $$

$$ \mathbf{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} $$

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) \to (\mathbf{e}_1', \mathbf{e}_2')$, or conversely its transpose $\mathbf{P} = \mathbf{R}^T$, which transforms the new system into the initial system and is called the passage matrix $\mathbf{P}: (\mathbf{e}_1', \mathbf{e}_2') \to (\mathbf{e}_1, \mathbf{e}_2)$. The rotation matrix and the passage matrix are orthogonal: $\mathbf{P}\mathbf{P}^T = \mathbf{P}^T\mathbf{P} = \mathbf{I}$ (identity matrix).

The transformation rule for a vector $\mathbf{u}$ under a change of coordinate system, that is, under a rotation of the basis vectors ($\mathbf{e}_1, \mathbf{e}_2$) $\to (\mathbf{e}_1', \mathbf{e}_2')$, is as follows, in terms of the passage matrix $\mathbf{P}$:

$$ \mathbf{u} \to \mathbf{u}' = \mathbf{P}^T \mathbf{u} \qquad (15) $$

$$ \mathbf{u}' \to \mathbf{u} = \mathbf{P} \mathbf{u}' $$

The quantity $\mathbf{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 (it is, rather, a pseudo-vector).

Finally, the transformation rule for a second rank tensor $\mathbf{T}$ under a rotation of the basis vectors ($\mathbf{e}_1, \mathbf{e}_2$) $\to (\mathbf{e}_1', \mathbf{e}_2')$ is as follows, in terms of the passage matrix $\mathbf{P}$:

$$ \mathbf{T} \to \mathbf{T}' = \mathbf{P}^T \mathbf{T} \mathbf{P} \to \mathbf{T}' \to \mathbf{T} = \mathbf{P} \mathbf{T}' \mathbf{P}^T \qquad (16) $$

$$ T_{ij} \to T'_{kl} = P_{ik} T_{ij} P_{jl} \quad T'_{kl} \to 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$); we insert the vector transformation rules $\mathbf{v} = \mathbf{P} \mathbf{v}'$ and $\nabla h = \mathbf{P} \nabla h'$; and we also insert the tensor transformation rule $\mathbf{K} = \mathbf{P} \mathbf{K}' \mathbf{P}^T$. This yields:

$$ \mathbf{P} \mathbf{v}' = - (\mathbf{P} \mathbf{K}' \mathbf{P}^T) \mathbf{P} \nabla h' = -\mathbf{P} \mathbf{K}' \nabla h' \quad (17) $$

$$ \mathbf{v}' = -\mathbf{K}' \nabla h' $$

which is just Darcy's law in the new system ($\mathbf{e}_i'$), with the transformed conductivity $\mathbf{K}'$ as expected.

Principal System, Diagonalization

The transformation rules serve 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:

$$ \det(\mathbf{K} - \lambda \mathbf{I}) = 0 \Rightarrow \lambda = \begin{vmatrix} \lambda_1 \\ \lambda_2 \\ \lambda_3 \end{vmatrix} \quad (18) $$

$$ \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} \Rightarrow [\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 of 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 principal basis, tensor $K_{ij}$ is represented by a diagonal matrix $\mathbf{K}^*$ with $K_{11}^* = \lambda_1, K_{22}^* = \lambda_2, K_{33}^* = \lambda_3$. The conductivity matrix $\mathbf{K}$ in any given coordinate system can always be expressed in terms of the diagonalized matrix $\mathbf{K}^*$ using the relation $\mathbf{K} = \mathbf{P}^* \mathbf{K}^* \mathbf{P}^T$, where $\mathbf{P}^*$ contains the $\mathbf{e}_i^*$'s as column vectors. More details in Appendix A.1.1 concerning the principal values and axes for the 2D case.

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

2.4.2. Directional Conductivity and Conductivity Ellipse

If $K_{ij}$ is a second rank tensor, symmetric and positive definite (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}$). These are functions of the direction of Darcy velocity $\mathbf{v}$ or of the gradient $\mathbf{J} = -\nabla h$, respectively. Direction angles can be used: polar angle $\theta$ in 2D; spherical angles ($\theta, \phi$) in 3D. It can 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 diameters 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} \text{ or : } v_i = K_{ij} J_j \text{ (with summation on repeated indices)} \quad (19) $$

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_{\text{FLUX}}$) is a scalar quantity obtained by “projecting” Darcy’s law on the direction of the “flux”, as follows:

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

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

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

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_{\text{FLUX}} = (\mathbf{n}_F^T \mathbf{K}^{-1} \mathbf{n}_F)^{-1} \quad \text{and} \quad K_{\text{GRAD}} = \mathbf{n}_G^T \mathbf{K} \mathbf{n}_G \quad (22) $$

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_{\text{FLUX}}(\theta_F)$ and $K_{\text{GRAD}}(\theta_G)$.

Now, it can be shown by inserting the diagonalization $\mathbf{K} = \mathbf{P}^*\mathbf{K}^*\mathbf{P}^T$ in these expressions that each of the quantities, $\sqrt{K_{\text{FLUX}}(\theta_F)}$ and $1/\sqrt{K_{\text{GRAD}}(\theta_G)}$, is the polar radius describing an ellipse whose principal axes are the eigenvectors of $\mathbf{K}$, and whose principal radii are $\sqrt{K_{ii}^*}$ (for the directional flux conductivity ellipse) and $1/\sqrt{K_{ii}^*}$ (for the directional gradient conductivity ellipse). Similarly, in 3D, $\sqrt{K_{\text{FLUX}}(\theta_F, \phi_F)}$ and $1/\sqrt{K_{\text{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. [40] designed 2D numerical flow experiments under imposed hydraulic gradient $\mathbf{J}_0$ in a synthetic fractured rock; they computed the directional $K_{\text{GRAD}}$ along the gradient direction, and then rotated $\mathbf{J}_0$ to obtain $K_{\text{GRAD}}$ for different angular directions, plotted the resulting $1/\sqrt{K_{\text{GRAD}}(\theta_G)}$, and fitted the resulting plot to an ellipse. Applications of the ellipse representation of conductivity anisotropy can also be found in forthcoming sections of the present paper; see the ellipses of equivalent block conductivity in Figures 15 and 16.

2.4.3. Non-Symmetric Permeability Tensor

Some of the upscaling techniques discussed in this paper lead to non-symmetric permeability tensors. This raises several questions about the validity of these methods. Before discussing this point further, let us first discuss the impact of non-symmetric permeability tensors on the Darcy velocity and derive some algebraic relations between anti-symmetry and diagonalization.

The degree of anti-symmetry can be separated algebraically from the other characteristics of the $\mathbf{K}$ tensor, such as its degree of anisotropy and the angles of its principal axes (e.g., directional permeability ellipse). Indeed, defining the symmetric part as $\mathbf{K}_S = (\mathbf{K} + \mathbf{K}^T)/2$, and the anti-symmetric part as $\mathbf{K}_A = (\mathbf{K} - \mathbf{K}^T)/2$, one can choose to diagonalize the symmetric part as usual, in an orthogonal basis (the resulting principal axes are orthogonal
and the resulting eigenvalues are real and positive). It can be shown that the rotation that diagonalizes $K_S$ has no effect on the anti-symmetric part $K_A$. Thus, the quadratic form $Q = x^T K x$ is seen to depend only on $K_S$:

$$Q = x^T K x = x^T (K_S + K_A) x = x^T K_S x + \frac{1}{2} x^T (K - K^T) x \quad (23a)$$

$$Q = x^T K_S x + \frac{1}{2} (x^T K x - x^T K^T x) = x^T K_S x + \frac{1}{2} (x^T K x - (x^T K^T x)^T) \quad (23b)$$ $$Q = x^T K_S x + \frac{1}{2} (x^T K x - x^T K^T x) = x^T K_S x$$

Similarly, it can be shown that $K$ is non-negative if $K_S$ is non-negative, while $K_A$ itself is always indefinite by construction. Therefore, the antisymmetric part $K_A$ cannot induce negativity of the $K$ tensor, no matter how anti-symmetric it might be. In conclusion, the symmetric part of the $K_{ij}$ tensor, $K_S = (K + K^T)/2$, can be analyzed in the same way as in the previous section (i.e., in terms of directional conductivities and ellipses), because the quadratic form “Q” of $K_S$ and of the non-symmetric $K$ are identical.

On the other hand, the anti-symmetric part, $K_A = (K - K^T)/2$, has a zero quadratic form. Therefore, the physical role of $K_A$ cannot be analyzed in terms of directional conductivities or conductivity ellipses. The full non-symmetric tensor $K = K_S + K_A$ can be analyzed more directly in 2D or 3D by diagonalizing it in two distinct non-orthogonal bases, that is, two distinct sets of principal axes and eigenvalues. The two distinct principal bases are related to the two distinct directions of the velocity and hydraulic gradient in the (non-symmetric) tensorial form of Darcy’s law. For details and interpretation, see for instance [69] (their Section 4.6 on Generalized principal axes). We propose below a physically based analysis of Darcy flow with non-symmetric tensorial conductivity.

The hydraulic effects of a non-symmetric equivalent permeability tensor $K_{ij}$ can be illustrated for the case of a 2D equivalent continuum as follows:

$$V = \begin{pmatrix} V_1 \\ V_2 \end{pmatrix} = - \begin{pmatrix} K_{11} \partial H / \partial x_1 + K_{12} \partial H / \partial x_2 \\ K_{21} \partial H / \partial x_1 + K_{22} \partial H / \partial x_2 \end{pmatrix} = \begin{pmatrix} K_{11} J_1 + K_{12} J_2 \\ K_{21} J_1 + K_{22} J_2 \end{pmatrix} \quad (24)$$

The above Darcy velocity vector is expressed for a given hydraulic gradient $J$, with $K_{21} \neq K_{12}$. Decomposing the velocity itself into symmetric and non-symmetric parts yields:

$$V = V_S + V_A \quad (25)$$

$$V_S = \begin{pmatrix} V_{1,S} \\ V_{2,S} \end{pmatrix} = \begin{pmatrix} K_{11} J_1 + 0.5(K_{12} + K_{21}) J_2 \\ 0.5(K_{12} + K_{21}) J_1 + K_{22} J_2 \end{pmatrix} = \begin{pmatrix} K_{11} J_1 + S J_2 \\ S J_1 + K_{22} J_2 \end{pmatrix} \quad (26)$$

$$V_A = \begin{pmatrix} V_{1,A} \\ V_{2,A} \end{pmatrix} = \begin{pmatrix} 0 \times J_1 + 0.5(K_{12} - K_{21}) J_2 \\ 0.5(K_{21} - K_{12}) J_1 + 0 \times J_2 \end{pmatrix} = \begin{pmatrix} +A J_2 \\ -A J_1 \end{pmatrix} \quad (27)$$

denoting $S = K_{12} + K_{21}$ and $A = K_{12} - K_{21}$.

Accordingly, the symmetric and antisymmetric contributions to velocity can be schematically represented as shown in Figure 3. In these schematics, the hydraulic gradient is considered as the forcing condition (input), and the Darcy velocity is the result (output). For illustration purposes, we consider the principal system ($x^*$) oriented at $\theta = \pi/4$ and the following non-symmetric $K = K_S + K_A$ tensor:

$$K_S = \begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix} = R^T K * R; K* = \begin{pmatrix} 3 & 0 \\ 0 & 1 \end{pmatrix} = R K_S R^T \quad (28)$$

$$R = \begin{bmatrix} \cos(\pi/4) & -\sin(\pi/4) \\ \sin(\pi/4) & \cos(\pi/4) \end{bmatrix}; R^T R = I \quad (29)$$

$$K_A = \begin{pmatrix} 0 & +1 \\ -1 & 0 \end{pmatrix} = \begin{pmatrix} 0 & +A \\ -A & 0 \end{pmatrix}; \Rightarrow K = K_S + K_A = \begin{pmatrix} 2 & 0 \\ -2 & 2 \end{pmatrix} \quad (30)$$
where:

* **R** is the passage matrix from reference frame to principal frame (and it is a $\pi/4$ rotation matrix);

* $\theta(x^* \to x) = +\pi/4$ is the angle of rotation from principal frame ($x^*$) to reference frame ($x$); and

* the reverse angle passing from reference frame to principal frame, ($x, x^*$), is $\theta(x \to x^*) = -\pi/4$.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a–f). (a,b): Schematic illustration of symmetric behavior, illustrating the symmetric contribution for a given hydraulic gradient **J** (orange arrow). The resulting Darcy velocity vector **V** is shown as a solid blue arrow. Here, the principal system is at an angle $\theta = -\pi/4$ with respect to the reference frame. (c,d): Schematic illustration of anti-symmetric behavior, illustrating the antisymmetric contribution for a given hydraulic gradient **J** (orange arrow). The Darcy velocity vector **V** is shown as a solid blue arrow. The orientation of the principal system ($\theta$) plays no role in this antisymmetric component of the flow. The antisymmetric deflection is always orthogonal (90 degrees). (e,f): Schematic illustration of the complete behavior of an anti-symmetric and inclined permeability tensor $\mathbf{K} = \mathbf{K}_S + \mathbf{K}_A$ for a given hydraulic gradient **J** (orange arrow, visible in (e), but hidden behind **V** in (f)). The resulting Darcy velocity vector **V** is shown as a solid blue arrow. The principal system of the symmetric part $\mathbf{K}_S$ is at an angle $\theta = -\pi/4$ with respect to the reference frame. The antisymmetric part ($\mathbf{K}_A$) exerts an orthogonal deflection (always at 90 degrees), which does not depend on the $\theta$ angle.</figcaption></figure>

In this example, the symmetric off-diagonal component is $S = [\mathbf{K}_S]_{12} = -1$, the antisymmetric part is $A = 2$, and the Frobenius norm is $||\mathbf{K}|| = ||\mathbf{K}_S|| = \sqrt{3 \times 3 + 1 \times 1} = \sqrt{10}$ (see Equation (134) for the Frobenius norm). Thus, in this particular example, the symmetric and antisymmetric parts are of the same order of magnitude (as occurs also in some of our experiments on equivalent block permeability).

Table 1 and Figure 3a–f show two examples illustrating the relative effects of non-diagonality and antisymmetry on the Darcy velocity vector **V** for a given hydraulic gradient vector **J**. The main result shown in these figures is that the antisymmetric part ($\mathbf{K}_A$) exerts a systematic orthogonal deflection (always at 90 degrees) that adds up to the effect of the symmetric part (compare Figure 3c,d).
<table><caption>**Table 1.** Antisymmetry effects on flow direction for a given hydraulic gradient vector **J**.</caption><thead><tr><th>Equations (Darcy's law)</th><th>Experiment No.1<br>Horizontal Gradient<br>[J<sub>1</sub>, J<sub>2</sub>] = [1, 0]</th><th>Experiment No.2<br>Vertical Gradient<br>[J<sub>1</sub>, J<sub>2</sub>] = [0, 1]</th></tr></thead><tbody><tr><td><b>Symmetric part:</b><br>V<sub>S</sub> = (2J<sub>1</sub> - J<sub>2</sub>)<br>(-J<sub>1</sub> + 2J<sub>2</sub>)</td><td>V<sub>S</sub> = (+2)<br>(-1)</td><td>V<sub>S</sub> = (-1)<br>(+2)</td></tr><tr><td><b>Antisymmetric part:</b><br>V<sub>A</sub> = (+J<sub>2</sub>)<br>(-J<sub>1</sub>)</td><td>V<sub>A</sub> = (0)<br>(-1)</td><td>V<sub>A</sub> = (+1)<br>(0)</td></tr><tr><td><b>Full:</b><br>V = (2J<sub>1</sub>)<br>(-2J<sub>1</sub> + 2J<sub>2</sub>)</td><td>V = (+2)<br>(-2)</td><td>V = (0)<br>(+2)</td></tr></tbody></table>

## 2.5. Boundary Conditions

Solving the equations of the local heterogeneous flow model and those of the upscaled macroscopic flow model requires in both cases to define a set of boundary conditions. In the following, we review the different conditions that are often considered and that we will use later.

### 2.5.1. Permeameter Boundary Conditions (PRM)

The permeameter boundary conditions correspond to the boundary conditions applied in the laboratory when measuring the hydraulic conductivity with a classical permeameter. They have been very often used to compute numerically the equivalent conductivity of heterogeneous media. The shape of the domain must include two parallel faces and correspond typically either to cylindrical cores or rectangular shaped samples. Based on this geometry, a constant head is prescribed on one of the two parallel faces and another constant head is prescribed on the opposite face, such that a global gradient is imposed in this direction. All the other faces are no-flow boundaries.

Note that, physically, laboratory samples are 3D (e.g., a cylindrical core, or a parallelepiped rectangle). However, quasi-2D samples can also be considered in the laboratory (thin slices, Hele-Shaw cells, etc.). To simplify our presentation, let us consider a two-dimensional experiment. We can express the permeameter conditions (PRM) as follows for a flow experiment with a head gradient imposed in the $x_1$ direction:

$$ 
\begin{array}{l}
h(0, x_2) = h_0, h(L_1, x_2) = h_0 + \Delta h_1, x_2 \in [0, L_2] \\
\mathbf{v}(x_1, 0) \bullet \mathbf{n} = \mathbf{v}(x_1, L_2) \bullet \mathbf{n} = 0, x_1 \in [0, L_1]
\end{array}
\qquad (31) $$

$L_1$ and $L_2$ represent the dimension of the edges of the block, $h_0$ the imposed head, $\Delta h_1$ the head difference between the two opposite faces, which are orthogonal to $x_1$. In this "$x_1$ flow experiment", a head gradient $\Delta h_1/L_1$ is imposed in the $x_1$ direction. The two other faces, orthogonal to $x_1$, are no-flow boundaries.

### 2.5.2. Dual Permeameter Conditions (PRM')

From a numerical point of view, it may be interesting to consider a novel configuration that generalizes the classical PRM permeametric conditions described above. Briefly, the new configuration consists of:

* two opposite faces with constant and identical non-zero fluxes, and

* linearly varying heads imposed on all other faces of the sample (instead of zero fluxes).

We propose to name these boundary conditions *dual permeameter* (PRM') boundary conditions.

### 2.5.3. Constant Gradient Conditions (GRD) from Linearly Varying Boundary Heads

A constant gradient condition (GRD) is imposed by prescribing linearly varying heads $h(x_1, x_2, x_3)$ along the boundary of the 2D domain or 3D domain, which may have complex
geometric shape (usually convex). The linear form $h(x_1, x_2, x_3)$ is defined in terms of a specified global or “far field” hydraulic gradient. This far field gradient occurs in the exterior of the porous domain and at its boundaries, hence the name “hydraulic immersion conditions” as if the porous sample were immersed in an ambient fluid characterized by a frozen macroscale gradient field. In the 3D case, for instance, this condition leads to a tri-linear macroscale head distribution on the boundary of the porous domain:

$$h(\mathbf{x}) = h_0 - \mathbf{J} \cdot \mathbf{x} \qquad (32)$$

where $\mathbf{J}$ [-] is sometimes called “hydraulic gradient”; it is the opposite of the 3D head gradient vector $\nabla h$ (macroscale), and it can be expressed as:

$$\mathbf{J} = \begin{bmatrix} J_1 \\ J_2 \\ J_3 \end{bmatrix} = - \begin{bmatrix} \partial h / \partial x_1 \\ \partial h / \partial x_2 \\ \partial h / \partial x_3 \end{bmatrix} \qquad (33)$$

With these notations, for any chosen orientation of the vector $\mathbf{J}$, the head distribution on the boundary is:

$$h(x_1, x_2, x_3) = h_0 - J_1 x_1 - J_2 x_2 - J_3 x_3 \qquad (34)$$

The same boundary conditions are applied on the microscale or macroscale experiments; therefore, we can use either the notation for micro $h(\mathbf{x})$ or macroscale $H(\mathbf{x})$ for the boundary conditions. In particular, for the case of a rectangular box-shaped domain, $J_1$ is the total head drop between any points at location $x_2 = x_3 = 0$ taken on the two opposite faces orthogonal ($\perp$) to axis $x_1$:

$$\mathbf{J} = \begin{bmatrix} J_1 \\ J_2 \\ J_3 \end{bmatrix} = \begin{bmatrix} \Delta h_1 / L_1 \\ \Delta h_2 / L_2 \\ \Delta h_3 / L_3 \end{bmatrix} \qquad (35)$$

The linear head distribution (Equation (34)) can be applied on any domain geometry in a numerical model while they are much more difficult to apply in the laboratory, as opposed to permeameter type boundary conditions. Linearly varying heads were used by many authors, including [40,43,78–81]. For instance, the last authors [81] applied linear head immersion conditions or equivalently constant gradient conditions (named here GRD) to a sample of fractured porous medium.

### 2.5.4. Piecewise Gradient Conditions (GRD')

As an alternative to the above GRD conditions, one may consider (and we will consider) piecewise constant gradient conditions (GRD') implemented as piecewise linear head conditions on the boundary. The GRD' conditions can be used to obtain analytically the local flow field on simple composite media (layered media, fractured media), and then perform averages to obtain explicitly an “exact” expression for the equivalent anisotropic conductivity tensor [80,82–84]. The last reference [84]), provides details in its Section 3.2.1 entitled “First step: equivalent permeability of unit blocks (upscaling single fracture/matrix blocks)”.

### 2.5.5. Constant Flux Immersion (FLX)

With permeameter boundary conditions (PRM), presented earlier, the flow is confined by impervious lateral boundaries around the sample. The flux immersion method (FLX) offers a more flexible manner to force the flow in any given direction, and this, for any domain shape. A constant macroscale vectorial flux $V_{BC}$ is imposed by prescribing
normal fluxes ($V_{BC} \bullet n$) on all elements of the boundary (on all points of the boundary $\Gamma$). Accordingly, let:

$$V_{BC} = (V_1 \ V_2 \ V_3)^T, \forall \mathbf{x} \in \Gamma : v_n(x) = V_{BC} \bullet n(x) \quad (36)$$

where $n(x)$ is the unit normal vector defined on any point $X$ of the boundary $\Gamma$ (outgoing normal). In addition, the head must be prescribed at a single point of the boundary to ensure unicity of the hydraulic head solution. These conditions can be understood as an immersion of the block in a frozen specific discharge vector field $V$ (macroscale). This vector is projected normally on each boundary point and used as a boundary condition. These boundary conditions have been used, for example, by Pouya and colleagues [43,85] for fracture networks.

The flux immersion condition, as described here, is a dual version of the GRD conditions discussed earlier. Note: by the same token, piecewise flux conditions (FLX') can also be constructed as an interesting analogy to the piecewise constant gradient conditions (GRD').

### 2.5.6. Periodic Conditions (PRD)

Periodic conditions have been used widely used in the framework of homogenization theory [12,37,74]. They are applied on regular domain geometry (rectangle, hexahedron, etc.). In the case of a 2D rectangular domain, we can express them as:

$$h(L_1, x_2) = h(0, x_2) - L_1 J_1 \quad (37)$$

$$h(x_1, L_2) = h(x_1, L_2) - L_2 J_2$$

$$v(x_1, 0) \bullet n = v(x_1, L_2) \bullet n$$

$$v(0, x_2) \bullet n = v(L_1, x_2) \bullet n$$

where $J_1$ and $J_2$ are the components of the vector $J$. In addition, a value of the hydraulic head must be imposed on a single boundary point, otherwise the problem would have an infinite number of solutions. Note that in numerical codes, e.g., with discretized finite volume formulations and iterative matrix solvers, this may not be needed, as the initial guess can serve as forcing head condition (this is the case for example with the BigFlow 3D code, documented in [86]).

Note that periodic boundary conditions on a rectangular domain can be interpreted as a formulation of the flow problem on a torus (doughnut-shaped domain). This interpretation also implies a periodic pattern of all material properties involved: the internal geometry of the medium is repeated periodically (or equivalently, it is mapped onto the surface of a torus). Thus, for a fracture network, this would require a periodic distribution of all conductive objects (fractures). Similarly, for a multilayered medium, this would require a periodic distribution of all layers.

### 2.5.7. Skin Immersion (SKN) or Border Region

It is possible to prescribe the boundary conditions not directly on the boundary of the domain of interest, but further away, by considering that the domain of interest is immersed within a larger domain having the same kind of heterogeneity. This method was proposed by [87], and used later by several authors such as [55].

### 2.5.8. General Case with Mixed Head and Flux Conditions

More generally, the boundary $\Gamma$ of the domain $\Omega$ may be divided into several parts, with various distributions of heads or fluxes on them. (Note: the term “flux” is sometimes used in this text to designate the specific discharge, or Darcy velocity vector.) For example, the boundary $\Gamma$ can be separated into two parts, $\Gamma = \Gamma_1 \cup \Gamma_2$, on which the head $h$ and flux $V$ are prescribed as given functions of space (respectively, $f$ and $g$):

$$\left| \begin{array}{l} \forall \mathbf{x} \in \Gamma_1 : h(\mathbf{x}) = f(\mathbf{x}) \\ \forall \mathbf{x} \in \Gamma_2 : v(\mathbf{x}) \bullet \mathbf{n} = g(\mathbf{x}) \end{array} \right. \quad (38)$$
### 3. Analytical Solutions of Darcian Flow Models

In this section, we present some analytical solutions of local-scale flow problems under specified boundary conditions in heterogeneous *continua* (these solutions are also applicable to macroscale problems with slowly evolving continuous permeability). In some cases, these analytical or quasi-analytical solutions have been backed up by numerical solutions obtained on fine grids (not detailed here). Such exact solutions can be useful for numerical benchmarks, and also for testing equivalent permeability criteria such as those studied in this paper. Note also that the equivalent permeabilities corresponding to other exact analytical solutions will be calculated explicitly for composite or layered media using our different equivalent permeability criteria (Section 5).

#### 3.1. Velocity-Based Vectorial Flow Equation with Heterogeneous Permeability

This section is a more detailed continuation of the velocity-based flow equations that were briefly introduced in Section 2.2 on local flow models (Section 2.2.2: velocity based flow equations). We start with an introduction concerning the origins, motivations, and possible interest of such velocity-based flow equations before presenting it briefly.

It has been observed in the literature on Darcy flow that, for steady flow in a heterogeneous multidimensional medium, it is possible to use a system of PDE's directly governing the Darcy velocity vector $v$ (or the flux density $q$) taking into account both Darcy's law and mass conservation without any explicit appearance of the pressure gradient or head gradient. This has several consequences:

* ✓ The detailed Darcy velocity field $v(x)$ resulting from the flux-based equations may be different from that obtained from the pressure-based equations. In theory they should be identical, but differences will arise due to the different approximations used either by perturbation methods or numerically in the flux-based and the pressure-based formulations; see for example Ababou ([68], Section 4.3 therein) and Akpoji et al. [72] concerning the flux vector variance in the case of a randomly heterogeneous spatially correlated medium (spectral perturbation solutions).

* ✓ Using the flux vector equations, new flux solutions may be obtained potentially for various types of heterogeneous media, deterministic or random, either by analytical methods or numerically; in the case of a finite domain, boundary conditions on the flux vector must be considered, since the pressure gradient does not appear explicitly in the flux equations (as will be seen).

* ✓ Another consequence of the flux-based approach is that it may lead to yet other definitions of the equivalent or effective permeability at the macroscale. Darcy's law must still be used, but the averaging is performed on the detailed flux or Darcy velocity field resulting from the flux-based equation (not from the pressure-based equation). Just for this reason, the resulting macroscale $\hat{K}_{ij}$ may be different.

Briefly, the “flux-based”, or Darcy velocity-based flow PDE, is shown below for a multi-dimensional porous medium with spatially variable hydraulic conductivity $K(x)$:

$$ \nabla^2 v + \nabla \times (\nabla LnK \times v) = 0 \qquad (39) $$

where the symbol $\times$ denotes the vector product.

Let us also denote $F(x)$ the log-permeability or log-conductivity field:

$$ F(x) = Ln\{K(x)\} \qquad (40) $$

The term in parentheses in Equation (39) can be interpreted as the vorticity vector $\Omega$ of the flow:

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

hence, the velocity-vorticity flow equation is:

$$ \nabla^2 v + \nabla \times \Omega = 0 \qquad (42) $$
A more advanced interpretation is developed in the aforementioned Appendices. For example, using indicial notations (with implicit summation on repeated indices), we have equivalently:

$$ \nabla^2 v_i = 2 \frac{\partial R_{ij}}{\partial x_j}; \quad R_{ij} = \frac{1}{2} \left( \frac{\partial v_i}{\partial x_j} - \frac{\partial v_j}{\partial x_i} \right) \qquad (43) $$

where $R_{ij}$ is the strain rate tensor. Note that this equation is a second order vector PDE (a system of 3 s order PDE's in 3D) governing the Darcy velocity vector (or "flux") in a heterogeneous porous medium. It is somewhat analogous to a vector Poisson equation or Helmholtz equation. In the present case, for Darcy flow, it expresses that the vector Laplacian of velocity is driven by vorticity or by the divergence of strain rate. This velocity PDE could be solved, without using the hydraulic head field, under conditions of a prescribed Darcy velocity $v$ along the boundary of the domain.

To sum up: the previous formulation follows the vectorial velocity-based formulation of Darcy flow formulated in Ababou [68] for a heterogeneous continuum in general. For details with particular emphasis on the case of random filed log-permeability, Ababou [68], Section 4.3 therein: "New Closed Form Perturbative Solution for the Flux Spectrum". This direct approach to velocity or flux solutions leads to stochastic analyses of flux variances, as in Akpoji, Ababou, De Smedt [70], based directly on the stochastic velocity equation. The velocity or flux-based equation was also presented in Zijl and Nawalany [69]; their Section 3.5 on "Directly calculated velocity components" (pp. 55-63) closely follows Ababou [68]'s work using different notations. We provide more details on the development of the velocity-based flow PDE in Appendix A.2: "Derivation of velocity-based PDE governing Darcy flow", which is based on Ababou [68] (Section 4.3 therein).

## 3.2. Philip 1986 Analytical Flow Solution for a Continuous Periodic K(x,y) Pattern

In his paper, Philip [88] examines a special case of saturated flow in 2D periodic heterogeneous porous medium such that the flow can be described quasi-analytically by simple, periodic flow functions: the potential function $\varphi(x, y)$ and the stream function $\psi(x, y)$. The potential $\varphi(x, y)$ is associated with the hydraulic head $h(x, y)$ and the stream function $\psi(x, y)$ is associated with the streamlines of the Darcy velocity field $v(x, y)$. He also discusses an extension in 3D, but here, we will briefly describe the 2D periodic flow pattern. In Philip's paper, the corresponding effective permeability is only mentioned in passing. Nevertheless, this particular flow field may be of interest for future equivalent permeability analysis.

Briefly, the 2D periodic permeability field $k(x, y)$ treated by Philip [88] is of the form:

$$ k(x,y)/k_0 = 1 + 2a \sin(x) \cos(y) + a^2 \cos^2(y) \qquad (44) $$

where $k_0$ is a reference value of permeability, obtained for instance at point $x = 0, y = \pi/2$ (and at all points located on horizontal lines $y = \pi/2$). In the remainder of this section, we take every where $k_0 = 1$, or equivalently, the permeability is rescaled by taking $k(x, y) \leftarrow k(x, y)/k_0$. With this scaling, the permeability is bounded below and above by: $0.5 \le k(x, y) \le 4.0$.

The corresponding Darcy velocity field $v(x, y)$ for this permeability pattern is given by:

$$ v_X = \partial\psi/\partial y = -k(x, y) \partial\varphi/\partial x; \quad v_Y = -\partial\psi/\partial x = -k(x, y) \partial\varphi/\partial y \qquad (45) $$

The stream function $\psi(x, y)$ was prescribed to be the following linear/periodic function:

$$ \psi(x, y) = y + a \sin(x) \cos(y); \quad 0 < a \le 1 \qquad (46) $$

Finally, the corresponding potential function $\varphi(x, y)$ can be determined by integration from:

$$ \partial\varphi/\partial y = a \cos(x) \sin(y) / \left\{ 1 + 2a \sin(x) \cos(y) + a^2 \cos^2(y) \right\} \qquad (47) $$
Note that $\partial\varphi/\partial y$ represents the local head gradient $\nabla h(x, y)$ in our notations. Finally, note that, in fact, the permeability pattern $k(x, y)$ of Equation (44) was inferred by Philip (1986) by solving exactly the following inverse problem: given the local flow pattern described by Equations (46) and (47), find the permeability field that satisfies this flow pattern: the result is the spatially distributed field $k(x, y)$ given just above in Equation (44).

Figure 4 displays the permeability pattern and the corresponding flow field of Philip (1986) for $a = 1$.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Illustration of the 2D periodic flow problem of Philip [70] for $a = 1$, in a square domain $[-\pi, +\pi] \times [-\pi, +\pi]$. **Top:** Permeability contours $k(x, y)$. **Bottom:** streamlines obtained as iso-values of the stream function $\psi(x, y)$.</figcaption></figure>
We now briefly consider an interpretation of the spatial pattern of $k(x,y)$ as a probability distribution (due to Philip 1986), and possible extensions of this analysis to obtain the effective macroscale permeability $K$ from knowledge of $k(x,y)$, $\varphi(x,y)$ and $[v_X(x,y), v_Y(x,y)]$.

i. Firstly, Philip [88] states that, for his periodic permeability pattern, the effective “apparent” permeability (which he denotes “K”) is by construction a scalar isotropic quantity equal to unity ($K^* = 1$). However, no details are given concerning the definition of such an apparent permeability.

ii. Secondly, Philip [88] was able to re-interpret the spatial distribution $k(x,y)$ of Equation (44) in terms of a probability density function of permeability (which he expresses in terms of an elliptic integral and plots in his Figure 3).

iii. We note that the (constant) effective permeability $K_{ij}$ (possibly tensorial) could be computed according to various definitions and criteria to be presented later in this paper. This could be done in two different ways:

a. by interpreting all local variables and coefficients in probabilistic terms (not only $k$, but also $v_X$, $v_Y$, and $\varphi$), and then use probabilistic averages $<...>$ to compute $K_{ij}$ such that:

$$ \langle v \rangle = -\mathbf{K} \langle \nabla \varphi \rangle \quad \text{or equivalently, with Einstein's indicial notation:} \\ \langle v_i \rangle = -\mathbf{K}_{ij} \langle \partial_i \varphi / \partial x_j \rangle \quad (i=1,2) $$

b. Alternatively, the same objective could be achieved by replacing the probabilistic averaging brackets above $<...>$ by spatial averaging brackets (being understood that the average is performed over a single periodic cell of the permeability pattern).

These suggestions have not been pushed further for the present periodic pattern; the probabilistic or spatial integrals needed are cumbersome, and it seems that the procedures suggested just above have not been carried out yet in the literature for this quasi-analytical flow pattern.

### 3.3. Analytical Velocity Field for Periodic Exponential-Cosine Permeability $K(x,y,z)$

The previous velocity-based equation (see also Appendix A.2) allows to obtain directly, in some cases, the Darcy velocity field for given permeability distributions, either exactly or approximately. We present below a velocity field in a finite domain with periodic permeability. This velocity is an approximation if one assumes a constant head gradient is imposed, but we will re-analyze it as an exact velocity solution with variable boundary head conditions on the finite size block. First, let us consider a 3D continuous porous block with a periodically variable permeability field (exponential-cosine function):

$$ K(x,y,z) = K_G \exp\{\sigma_{LnK} \cos(w_X x + w_Y y + w_Z z)\} \quad (48) $$

where $(w_X, w_Y, w_Z)$ is the spatial frequency or wave vector, $\sigma_{LnK}$ is the amplitude of Ln(K) fluctuations, and $K_G$ is the geometric mean permeability. The following analytical velocity field satisfies exactly 3D steady state mass conservation $\mathrm{div}(\vec{v}) = 0$, and it also satisfies Darcy's law—albeit approximately (the approximation is first order in $\sigma_{LnK}$ and it improves as $\sigma_{LnK}$ decreases and becomes much less than unity).

$$ v_X(x,y,z) = \bar{V}_X \left\{ 1 + \sigma_{LnK} \left( 1 - \frac{w_X^2}{w_X^2 + w_Y^2 + w_Z^2} \right) \cos(w_X x + w_Y y + w_Z z) \right\} \quad (49) $$

$$ v_Y(x,y,z) = -\bar{V}_X \sigma_{LnK} \left( \frac{w_X w_Y}{w_X^2 + w_Y^2 + w_Z^2} \right) \cos(w_X x + w_Y y + w_Z z) $$

$$ v_Y(x,y,z) = -\bar{V}_X \sigma_{LnK} \left( \frac{w_X w_Y}{w_X^2 + w_Y^2 + w_Z^2} \right) \cos(w_X x + w_Y y + w_Z z) $$
**Remark 1.** The above solution corresponds to a periodic medium with a single wave vector. It was used by one of us (R.A.) circa 1992 to test Lagrangian Particle Tracking algorithms on massively parallel Connection Machine. A more general Fourier expansion approximation based on spectral perturbation theory [32,68] was also developed for a broader spectrum of wave-vectors, and a similar general approximation was independently presented and used by Kapoor and Kitanidis [89] for illustrating solute dilution and concentration fluctuations in a heterogeneous aquifer (cf. their appendix: "Model Periodic Hydraulic Conductivity and Flow Field", p. 313 therein).

Figure 5a shows a plot of the 2D periodic exponential-cosine permeability field:

$$K(x, y) = K_G \exp\{\sigma_{LnK} \cos(w_X x + w_Y y)\} \quad (50)$$

with $w_X = \frac{2\pi}{\lambda_X}$; $w_Y = \frac{2\pi}{\lambda_Y}$; $w_Z = 0$; $\lambda_X = \lambda_Y = 1$; $K_G = 1$; $\sigma_{LnK} = Ln10 \approx 2.30$.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a,b) Analytic Darcy flow problem with exponential-cosine permeability: permeability field, streamlines, and hydraulic head field, from top to bottom. (a): contours of the periodic permeability field $K(x, y)$ with $\sigma_{LnK} \approx 2.30$; the bold black contours correspond to geometric mean permeability $K_G = 1$; the broad white zones without contour lines correspond to the relatively flat minima of $K(x, y)$. (b): streamlines of the analytical Darcy velocity field ($v_X(x, y)$, $v_Y(x, y)$) given in the text.</figcaption></figure>

Figure 5b shows a plot of this 2D periodic permeability field, and of the streamlines corresponding to the 2D analytical Darcy velocity field ($v_X(x, y)$, $v_Y(x, y)$) with the above parameters and $(\bar{V}_X, \bar{V}_Y) = (1, 0)$.

Before proceeding further, let us note that the previous velocity field solution was inspired by spectral perturbation solutions of Darcy flow in an infinite domain with randomly heterogeneous permeability (the periodic permeability proposed here represents only a single non-random Fourier mode from the spectral perturbation approach). The resulting periodic velocity field above (Equation (49)) is only an approximation of the local
velocity if one assumes that the flow domain is infinite and that the flow is submitted to a constant large scale hydraulic gradient: in that case, the velocity of Equation (49) satisfies exactly mass conservation but only approximately Darcy's law.

However, let us now take a different viewpoint in two ways: (i) by restricting the flow problem to a finite domain (e.g., 2D square or rectangular domain), and (ii) by seeking boundary conditions such that the local velocity field given above of Equation (49) becomes exact. This is achieved by the following “inverse” procedure, that is, inverting Darcy’s law as follows:

$$j(x, y) = v(x, y)/K(x, y) \qquad (51)$$

where $j(x, y)$ is the local hydraulic gradient defined by $j(x, y) = -\nabla h(x, y)$, and $K(x, y)$ is the periodic “exponential-cosine” function of space, defined earlier.

Note that the flow field defined by $v(x, y)$, $j(x, y)$, and $h(x, y)$ now satisfies exactly both mass conservation ($div(v) = 0$) and Darcy's law ($v = -K(x, y)\nabla h$) with the specified exponential-cosine permeability field $K(x, y)$.

It is now possible to explicitly obtain the head distribution $h(x, y)$ inside the square domain, and along its boundary as well (four faces: left, right, bottom, top). To obtain $h(x, y)$ at any interior or boundary point $(x, y)$, one may first integrate horizontally the horizontal gradient component $j_X(x, y)$ and then vertically the vertical gradient component $j_Y(x, y)$, following for instance the integration path illustrated in the schematic of Figure 6.

<figure><img src="image_9.png" /><figcaption>Figure 6. Schematic representation of an integration path to obtain the head $h(x, y)$ at any point $(x, y)$, given the hydraulic gradient field $j(x, y) = -\nabla h(x, y)$.</figcaption></figure>

Assuming for convenience a rectangular domain $[0, L_X] \times [0, L_Y]$, the quasi-analytical result is:

$$h(x, y) = h(0, 0) - \int_{Y=0}^{Y=y} j_Y(0, Y) dY - \int_{X=0}^{X=x} j_X(X, y) dX \qquad (52)$$

where the vector $j(x, y)$ is known analytically from the previous inverse Darcy relation. The simple integrals were carried out with a simple histogram method (mimicking Rieman integral), and the resulting head distribution inside the domain is displayed in Figure 7 for the same case of highly variable periodic $K(x, y)$ with $\sigma_{LnK} = Ln10 \approx 2.30$ as in the previous figures.
**Remark 2.** $\sigma_{LnK}$, the standard deviation of LnK(x,y), plays the role of an amplitude parameter here; the heterogeneous domain is very variable for $\sigma_{LnK} \gg 1$ and, on the contrary, it becomes nearly homogeneous for $\sigma_{LnK} \ll 1$.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Quasi-analytical hydraulic head field $h(x, y)$ given in the text, shown as bold contours. The contours of periodic permeability $K(x, y)$ with $\sigma_{LnK} = Ln10 \approx 2.30$ are also superimposed as thin curves. The broad white bands of $K(x, y)$ correspond to zones of minimal permeability, coinciding also with zones of large head gradient.</figcaption></figure>

Finally, from the previous head solution $h(x, y)$, which is valid both at interior points and at boundary points, we now deduce the four head profiles along the four boundary faces of the domain, as follows Equations (53)–(56):

$$ \text{Bottom face } y = 0: h(x, 0) = h(0, 0) - \int_{X=0}^{X=x} j_X(X, 0) dX \quad (53) $$

$$ \text{Top face } y = L_Y: h(x, L_Y) = h(0, 0) - \int_{Y=0}^{Y=L_Y} j_Y(0, Y) dY - \int_{X=0}^{X=x} j_X(X, L_Y) dX \quad (54) $$

$$ \text{Left face } x = 0: h(0, y) = h(0, 0) - \int_{Y=0}^{Y=y} j_Y(0, Y) dY \quad (55) $$

$$ \text{Right face } x = L_X: h(L_X, y) = h(0, 0) - \int_{Y=0}^{Y=y} j_Y(0, Y) dY - \int_{X=0}^{X=L_X} j_X(X, y) dX \quad (56) $$

The corresponding profiles of boundary heads are plotted in various ways, for the high variability case ($\sigma_{LnK} = Ln10 \approx 2.30$), in Figure 8. We have taken the head value at the lower left corner to be zero: $h(0, 0) = 0$. *Reminder*: in all such figures (velocities, streamlines, heads), the domain is a square with coordinates $x \in [-1, +1]$ and $y \in [-1, +1]$, and therefore, the domain size $L_X \times L_Y$ is $2 \times 2$.

Note in particular that the continuity of hydraulic head at the four corner nodes is satisfied; for instance, at the bottom right corner ($L_X, 0$), the head value from the bottom horizontal profile $h(x, 0)$ at $x = L_X$ is the same as the head value from the right vertical profile $h(L_X, y)$ at $y = 0$. They can be seen to be identical graphically from the plotted boundary profiles of Figure 8. We have verified indeed that head continuity is ensured
with relative precision better than 1/1000, and it is also continuous analytically. The unique value of head at the bottom right corner is numerically $h(L_X, 0) \approx -4.832$, and analytically:

$$h(L_X, 0) = - \int_{X=0}^{X=L_X} j_X(X, y) dX \quad (taking h(0,0) = 0).$$

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Plot of the 4 profiles of hydraulic head $h(x)$ or $h(y)$ along the 4 boundary faces of the square domain with periodic "exponential-cosine" permeability $K(x, y)$ for $\sigma_{LnK} = Ln10 \approx 2.30$ (high variability). It can be seen that the continuity of hydraulic head at the 4 corners of the square domain is ensured (more on this in the text).</figcaption></figure>

Furthermore, due to the periodic permeability field, the head field has a number of symmetries, as can be seen from $h(x, y)$ in Figure 7, and from the boundary head profiles of Figure 8. Thus, in Figure 8, the bottom profile $h(x)$ in dark blue is nearly identical to the left profile $h(y)$ in cyan color. Similarly, the top profile $h(x)$ in magenta is nearly identical to the right profile $h(y)$ in green. The slight discrepancies between the bottom/left profiles and between the top/right profiles are only due to coarse numerical integration (the graphics shown here use a discretization of the domain into 50 × 50 pixels or 51 × 51 nodes; these discrepancies will vanish for finer numerical integration grids).

In summary, an exact quasi-analytical finite domain flow under prescribed spatially variable head conditions was presented in this subsection for a periodic permeability field. The sample appears to be stratified with an angle of stratification of 45 degrees (in our example). This sample flow problem could be further exploited to infer the equivalent homogeneous permeability using one or several of the methods or criteria that will be presented later in this paper.

## 4. Equivalent Conductivity: The Criteria

**Definition 1.** The equivalent block conductivity $K$ is the hydraulic conductivity of a fictitious upscaled homogeneous medium such that certain average quantities (the equivalence criteria) are equal if we compute them both on the real heterogeneous medium and on the upscaled medium when the same boundary conditions are imposed on the two media. NB: more precisely, in some methods it is the equivalent conductivity matrix that is obtained; in other methods the resulting conductivity $K$ can be shown to be a true second rank tensor $K_{ij}$.
To apply the previous definition, the equivalence criteria must be defined precisely. In the next paragraphs, we describe the average quantities that have been used to define the equivalence between the heterogeneous and upscaled media. In each case, we distinguish (when necessary) between the three types of heterogeneous media defined earlier: (1) continuous, (2) composite, and (3) fracture networks. In the case of fracture networks, hydraulic conductivities vanish in most of the domain except in the fracture set, and for this reason, the permeable domain is partially disconnected, and so is its boundary Γ. Due to this extreme type of discontinuity, we use special dedicated formulations for flux averages and equivalent conductivities. In all cases, the previous formal definition requires the computation of the equivalence criteria both on the heterogeneous and on the homogeneous media.

## 4.1. Average Quantities and Equivalence Criteria

In this section, we present several criteria (for instance based on different kinds of averages of flux and hydraulic gradient, or total dissipated power, etc.) that can be used to define the equivalent block-scale permeability. For each type of criterion, the case of a porous continuum sample is treated in detail, and in most cases, the same expressions hold also for a composite porous sample (in addition, the corresponding expressions for the special case of a block consisting of a 2D fracture network are summarized).

### 4.1.1. Net Surface Flux (NSF)

The Net Surface Flux (NSF) criterion is the first and most obvious criterion that was used [90,91]. It is based on the analogy with the permeameter experiment.

The idea is that the net fluxes through the surfaces of the heterogeneous and homogeneous media must be identical. It is perfect for isotropic medium. However, it leads to some difficulties for anisotropic medium; when all the faces of the domain are open to flow, the flux entering the domain through a given face can be different from the flux leaving the domain through the opposite face. A way to circumvent this issue is to use an average: the mean between the fluxes on the two opposite faces (in the case of rectangular geometry).

#### NSF Criterion for a 2D Rectangular Porous Sample

Thus, for a 2D rectangular porous sample, the $x_1$ component of the net surface flux vector is then defined by the following expression for continuous and composite porous media:

$$ \langle \mathbf{v} \rangle_1^{\text{NSF}} = \frac{\int_{\Gamma_{B1}} \mathbf{v} \cdot \mathbf{n}_{\Gamma_{B1}} d\gamma - \int_{\Gamma_{A1}} \mathbf{v} \cdot \mathbf{n}_{\Gamma_{A1}} d\gamma}{2 L_2} \qquad (57) $$

where the minus sign results from the opposite orientation of the normal vectors ($\mathbf{n}_{\Gamma_{A1}} = -\mathbf{n}_{\Gamma_{B1}}$) for two opposite faces in rectangular geometry (in this paper, the normals are outgoing, see Figure 1).

#### NSF Criterion for a 2D or 3D Porous Sample of Arbitrary Shape

For the sake of generality, we propose a set of expressions that can be used on any geometry of the porous sample in two or three dimensions. For that purpose, we assume that only surface fluxes at the boundary of the domain are available to measurements and not vectors. Under these constraints, we derive (see Appendix A.4) the following set of equations in two dimensions.

$$ \langle \mathbf{v} \rangle_1^{\text{NSF}} = \frac{N_2 F_1 - N_{12} F_2}{N_1 N_2 - N_{12} N_{12}}, \quad \langle \mathbf{v} \rangle_2^{\text{NSF}} = \frac{N_1 F_2 - N_{12} F_1}{N_1 N_2 - N_{12} N_{12}} \qquad (58) $$

where $N_1$, $N_2$, and $N_{12}$ are three scalars that depend only on the geometry of the domain and on the coordinate system, and $F_1$ and $F_2$ are two oriented and integrated surface fluxes, as follows:

$$ \langle \mathbf{v} \rangle_1^{\text{NSF}} = \frac{N_2 F_1 - N_{12} F_2}{N_1 N_2 - N_{12} N_{12}}, \quad \langle \mathbf{v} \rangle_2^{\text{NSF}} = \frac{N_1 F_2 - N_{12} F_1}{N_1 N_2 - N_{12} N_{12}} \qquad (59) $$
$$N_{12} = \int_{\Gamma} (\mathbf{u}_1 \cdot \mathbf{n}) (\mathbf{u}_2 \cdot \mathbf{n}) d\gamma \qquad (60)$$

$$F_1 = \int_{\Gamma} (\mathbf{v} \cdot \mathbf{n}) (\mathbf{u}_1 \cdot \mathbf{n}) d\gamma, \quad F_2 = \int_{\Gamma} (\mathbf{v} \cdot \mathbf{n}) (\mathbf{u}_2 \cdot \mathbf{n}) d\gamma \qquad (61)$$

In these equations, $\mathbf{u}_i$ is the unit vector along direction $i$. In the special case of a rectangular domain, $F_i$ boils down to the total flux normal to the face orthogonal to $\mathbf{u}_i$. This set of equations leads to Equation (57) when the domain has a rectangular shape. The formula is consistent, that is, when applied to a constant velocity field within the domain it yields the components of this constant velocity vector (as it should). The general expressions for the 3D case are given in Appendix A.5.

### NSF Criterion for a 2D Fracture Network

Although we have not systematically specialized our analyses to fracture networks in this paper, it may be useful to consider the previous “NSF” criterion when applied to a 2D fracture network consisting of an impervious block traversed by thin fractures (segments in 2D). For such a fracture network, NSF expressions similar to the previous ones can be obtained via a discrete summation of the fluxes at the boundary of the fractured block.

To illustrate the principle of this summation, Figure 9 shows a simple square block. On the external boundary of the block, only the incoming or outcoming fluxes are assumed to be measurable. These quantities are defined for each node $i$ located on the boundary and denoted $q_i$. By convention, $q_i$ is positive if the flux enters the block and negative if the flux exists the block. With this notation, the horizontal component of the Net Surface Flux can be expressed as the average flux between the two opposite faces (see Appendix A.4 for a generalization).

$$\langle \mathbf{v} \rangle_{1}^{\text{NSF}} = \frac{q_{2} - q_{3}}{2L}, \quad \langle \mathbf{v} \rangle_{2}^{\text{NSF}} = \frac{-q_{4}}{2L} \qquad (62)$$

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Schematic illustration of the NSF criterion on a simple fracture network. Note that only the normal fluxes perpendicular to the faces of the block are available at the nodes that are located on the boundary (compare to the VSF method of Figure 10).</figcaption></figure>

Using either Darcy's law or Poiseuille's law, it is possible to express the vectorial fluxes $\mathbf{q}_B$ at all boundary nodes of the fracture network. One can then express formally the Net Surface Flux by multiplying the fluxes $\mathbf{q}_B$ by a matrix $\mathbf{P}$ having two lines (for two vector components in 2D) and as many columns as the number of boundary nodes. The coefficients of this matrix are computed in order to project the fluxes on the right direction, and we obtain the mean NSF velocity as follows in terms of border fluxes $\mathbf{q}_B$:

$$\langle \mathbf{v} \rangle_{\text{NSF}} = \mathbf{P} \mathbf{q}_B \qquad (63)$$
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Schematic illustration of the VSF criterion on a simple fracture network. Note that the directions of the fluxes are also available (contrary to the NSF method of Figure 9). In this case, the vectors are summed over all the edges touching the border.</figcaption></figure>

Conclusions on the NSF Criterion

Even if the Net Surface Flux (NSF) criterion has some interesting features, such as its consistency, it does not allow capturing properly the deviation of flow that can occur in a medium. For example, if the flux enters a face at the base of the face and leaves the opposite face near a corner, it is clear that there exists a diagonal component in the flux that we cannot capture with the NSF. This is why other formulas are required; this is the motivation for the formulas described in the next sections.

4.1.2. Vectorial Surface Flux (VSF)

In some situations, it is possible to directly access the vectorial fluxes at the boundary of the domain. For example, the flux vectors (or velocity vectors) can be accessed quite directly at the boundary of a fractured block consisting of a 2D fracture network, but more generally, they can also be computed on the boundary of any continuous or composite porous block.

VSF Criterion for a Porous Sample in 2D or 3D

It is quite natural to define an average vectorial surface flux at the boundary of the sample by simply averaging all those vectors. This can be done by a surface integral as follows, for a sample of continuous or composite porous medium (generally in three-dimensional space).

$$ \langle \mathbf{v} \rangle_{\text{VSF}} = \frac{1}{\Gamma} \int_{\Gamma} \mathbf{v} \, d\gamma \qquad (64) $$

VSF Criterion for a 2D Fracture Network

In the case of a 2D fracture network, one can use instead a discrete summation over all the individual fractures hitting the boundary of the block, as shown for illustration in the example of Figure 10:

$$ \langle \mathbf{v} \rangle_{\text{VSF}} = \frac{1}{2L} (\mathbf{v}_1 + \mathbf{v}_2 + \mathbf{v}_3 + 2\mathbf{v}_4) \qquad (65) $$
### 4.1.3. Volume Averaged Flux (VAF)

Another interesting alternative technique to get a mean vectorial flux in the block is to average the Darcy velocities over the whole domain of interest $\Omega$ as proposed by [92]. This is what we define as the Volume Average Flux

$$ \langle \mathbf{v} \rangle_{\text{VAF}} = \frac{1}{\Omega} \int_{\Omega} \mathbf{v} \, d\omega \qquad (66) $$

This concept can be less intuitive than the two previous averaging methods (NSF and VSF) when one thinks about a laboratory experiment in which the Darcy velocities inside the block are not easily accessible. Of course, this is not an issue for numerical experimentation. However, this is also not a problem in a laboratory experiment because the volume averages can be replaced by surface integrals [93] as shown in Appendix A.7:

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

Developing the previous equation on square geometries leads to an average of the total surface fluxes weighted by the coordinates of the centers of gravity of those fluxes on the different faces (Appendix A.7.3 and Table 2). Note that these equations have been used to construct a prototype of tensorial permeameter in the laboratory [25].

<table><caption>**Table 2.** Equivalent permeability criteria: specialization of volume averages to surface averages for specific boundary conditions and block geometries (see Appendix A.7 for proofs).</caption><thead><tr><th></th><th>Volume Average Gradient (VAG)<br>∇h<sub>VAG</sub></th><th>Volume Average Flux (VAF)<br>v<sub>VAG</sub></th><th>Total Dissipated Power (TDP)<br>−Pf/ρg</th></tr></thead><tbody><tr><td>Definition (via volume average)</td><td>1/Ω ∫Ω ∇h dω</td><td>1/Ω ∫Ω v dω</td><td>∫Ω (v•∇h) dω</td></tr><tr><td>Surface integral (inferred from the volume average)</td><td>1/Ω ∫Γ h n dγ</td><td>1/Ω ∫Γ x (v•n) dγ</td><td>∫Γ h v•n dγ</td></tr><tr><td>Gradient conditions GRD: constant gradient J. Any geometry</td><td>−J</td><td>1/Ω ∫Γ x (v•n) dγ</td><td>J• &lt; v &gt;<sub>VAF</sub></td></tr><tr><td>Expression for polyhedral geometry (for any BC)</td><td>1/Ω ∑<sup>n</sup><sub>i=1</sub> h<sub>i</sub>n<sub>i</sub>S<sub>i</sub></td><td>1/Ω ∫Γ x v•ndγ</td><td>∫Γ h v•n dγ</td></tr><tr><td>Expression for rectangular geometry (for any BC)</td><td>⟨∇h⟩<sub>i</sub> =<br>1/L<sub>i</sub>Γ<sub>Bi</sub> ∫Γ<sub>Bi</sub> h dγ<br>− 1/L<sub>i</sub>Γ<sub>Ai</sub> ∫Γ<sub>Ai</sub> h dγ</td><td>⟨v⟩<sub>i</sub> =<br>1/Ω (q<sub>Ai</sub>x<sub>i</sub><sup>Ai</sup> + q<sub>Bi</sub>x<sub>i</sub><sup>Bi</sup>)<br>+ q<sub>Ai</sub>x<sub>Gi</sub><sup>Ai</sup> + q<sub>Bi</sub>x<sub>Gi</sub><sup>Bi</sup>)</td><td>∫Γ h v•n dγ</td></tr><tr><td>Periodic conditions PRD</td><td>−J</td><td>⟨v⟩<sub>1</sub> = q<sub>B1</sub>/Γ<sub>1</sub>; ⟨v⟩<sub>2</sub> = q<sub>B2</sub>/Γ<sub>2</sub></td><td>J<sub>1</sub> q<sub>B1</sub> + J<sub>2</sub> q<sub>B2</sub></td></tr><tr><td>Permeameter conditions PRM (here for flow parallel to axis x<sub>1</sub>).<br>Rectangular geometry (2D example)</td><td>⟨∇h⟩<sub>1</sub> = Δh<sub>1</sub>/L<sub>1</sub> =<br>Δh<sub>B1</sub>−Δh<sub>A1</sub>/L<sub>1</sub><br>⟨∇h⟩<sub>2</sub> = 1/L<sub>2</sub>Γ<sub>B2</sub> ∫Γ<sub>B2</sub> h dγ<br>− 1/L<sub>2</sub>Γ<sub>A2</sub> ∫Γ<sub>A2</sub> h dγ</td><td>⟨v⟩<sub>1</sub> = q<sub>B1</sub>/Γ<sub>1</sub><br>⟨v⟩<sub>2</sub> = q<sub>B1</sub>(x<sub>G2</sub><sup>B1</sup>−x<sub>G2</sub><sup>A1</sup>)/Γ<sub>1</sub></td><td>Δh<sub>1</sub> q<sub>B1</sub></td></tr><tr><td>Flux immersion (FLX) with fixed constant flux V</td><td>1/Ω ∫Γ h n dγ</td><td>V</td><td>V•⟨∇h⟩<sub>VAG</sub></td></tr></tbody></table>

NB: some methods not shown in this table are VSF and NSF (directly involving surface fluxes).
## VAF Criterion for a 2D Fracture Network

On the example of the previous fracture network (see also Figure 2), applying the VAF criterion in a discrete manner yields (for a square domain):

$$ \langle \mathbf{v} \rangle_{\text{VAF}} = \frac{1}{L^2} (q_2\mathbf{x}_2 + q_3\mathbf{x}_3 + q_4\mathbf{x}_4) \qquad (68) $$

### 4.1.4. Volume Averaged Gradient (VAG)

Up to this point, we only discussed average quantities related to the Darcy velocities. However, depending on the boundary conditions, it may be necessary as well to estimate an average hydraulic gradient.

For example, when imposing permeameter type boundary conditions, only the head gradient in one direction is imposed during a single experiment, but the perpendicular gradient components also need to be computed.

As another example, when imposing a fixed constant vectorial flux to the boundaries of the domain, these fixed velocities are interpreted as the block scale velocities, and only the average head gradient requires to be computed as a function of the heterogeneity.

In all those situations, the average gradient is computed as the volume integral over the local gradients within the domain, as proposed by Rubin and Gómez-Hernández [92].

$$ \langle \nabla h \rangle_{\text{VAG}} = \frac{1}{\Omega} \int_{\Omega} \nabla h \, d\omega \qquad (69) $$

Like for the Volume Average Flux, one can transform the volume integral into a surface integral [93] to facilitate its computation or its measurement in a real experiment (see Table 2 and Appendix A.7.2).

Among the results concerning this quantity, it is important to note that the volume average gradient is equal to the far field gradient J imposed through linearly varying head boundary conditions for any geometry, and it is also true for periodic boundary conditions. This result is true both for the heterogeneous and homogeneous media because it depends only on the head distribution along the boundary.

In conclusion, the Volume Average Gradient (VAG) is identical to the opposite of the macroscale hydraulic gradient vector, which serves to impose linearly varying head conditions on the boundary of the domain. That is, we have:

$$ \langle \nabla h \rangle_{\text{VAG}} = -\mathbf{J} \qquad (70) $$

where J is the hydraulic gradient vector in the linearly varying head boundary conditions:

$$ h(\mathbf{x}) = h_0 - \mathbf{J} \cdot \mathbf{x} \qquad (71) $$

A similar linear distribution of hydraulic heads can be expressed at the boundary of a fracture network, by prescribing $h(\mathbf{x})$ at all fracture boundary points. These are the intersection points of fractures with the immaterial boundary of the computational domain (e.g., a rectangular block).

### 4.1.5. Total Dissipated Power (TDP)

The last very important physical quantity that we will consider in this paper is the Total Dissipated Power. This quantity was frequently used to define the equivalent hydraulic conductivity [30,94–96].

When flow occurs through a porous or fractured medium, a dissipation of energy occurs, which corresponds to the hydraulic head loss along the flow field. This energy loss is equivalent to the work of the viscous forces that resist the flow.
## TDP for a Porous Continuum

In the case of a porous continuum, the total energy dissipated per unit time within the porous domain is the total dissipated power P [Watt] defined as

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

where g is the acceleration of gravity [m/s²] and ρ is the mass density [kg/m³] of the fluid.

## TDP for a Fracture Network

Briefly, for a 2D fracture network, the equivalent expression of the Total Dissipated Power P is of the form:

$$P = \rho g \sum (-\mathbf{v} \cdot \nabla h) \Delta\omega \qquad (73)$$

where the summation is over all conductive arcs, or fracture segments, in the network.

## 4.2. Defining the Equivalent Conductivity Tensors

In this section, we describe how the equivalent block conductivity tensor (or at least the equivalent conductivity matrix) can be expressed and calculated, depending on the choice of the equivalence criteria and of the boundary conditions imposed on the sample block. We treat generally the case of continuous or possibly composite porous samples, and we also briefly provide the corresponding expressions for the special case of a 2D fractured block consisting of a network of fracture segments. We end up with a special consideration of the dual matching procedure for determining the equivalent permeability and resistivity, which, although not implemented numerically in this paper for the various block scale samples of Section 6, has proven useful in the context of upscaling for randomly heterogeneous media via analytical perturbation methods (e.g., [75]).

### 4.2.1. Equivalent K from Permeametric Experiment: Diagonal K Matrix (DIAG)

Since the first numerical studies of the equivalent permeability of heterogeneous media [91], the equivalent permeability “K” is often estimated by imposing permeameter type boundary conditions and computing the net surface flux through the heterogeneous medium in the direction of the imposed gradient. Assuming that the main axes of anisotropy of the medium are parallel to the axes of the sample, the net surface flux through the homogeneous medium is simply equal to the imposed head gradient times the diagonal element of the equivalent conductivity tensor corresponding to the direction of the flow experiment.

Consequently, one can estimate the equivalent (diagonal) conductivity matrix in the coordinate system of the porous sample, provided three numerical simulations (in 3D) or two numerical simulations (in 2D). Taking the 2D example yields the following diagonal K matrix (note: it is not the equivalent conductivity “tensor” that is “diagonal”; rather, it is the equivalent conductivity *matrix* K that is *diagonal* in the coordinate system of the block.):

$$K_{DIAG} = - \begin{pmatrix} \frac{\langle v_1 \rangle_{NSF}^1}{\langle \partial_1 h \rangle_{VAG}^1} & 0 \\ 0 & \frac{\langle v_2 \rangle_{NSF}^2}{\langle \partial_2 h \rangle_{VAG}^2} \end{pmatrix} \qquad (74)$$

In that expression, $\langle v_1 \rangle_{NSF}^1$ is the net surface flux in the $x_1$ direction computed from the first numerical experiment (hence the superscript 1) with a gradient $\langle \partial_1 h \rangle_{VAG}^1$ imposed in the $x_1$ direction. Similarly, $\langle v_2 \rangle_{NSF}^2$ is computed from the numerical experiment in the perpendicular direction.

The same expression can be written in 3D; it only requires running an additional flow experiment in the $x_3$ direction.

In the following, the DIAG method correspond to the results obtained with Equation (74) whether they are computed from experiments made using permeameter type or linearly varying head boundary conditions. In other words, the DIAG method consists of using the
NSF criterion and neglecting the off-diagonal components of the equivalent permeability matrix expressed in the coordinate system of the porous sample. This method is obviously not adequate for identifying the full permeability tensor in the general anisotropic case, but is used here to allow a comparison with the other techniques that are presented below.

Note: "matrix" and "tensor" should be distinguished. Indeed, assuming that the equivalent K is a second rank tensor ($K_{ij}$), then the permeability matrix K (diagonal or not) is merely the representation of this tensor in the reference system (in the present case, for permeametric conditions, the reference system is aligned with the porous block sample).

### 4.2.2. K from Net Surface Flux (NSF)

A natural extension of the previous definition (DIAG) to anisotropic media is to consider the Net Surface Flux (NSF) to define an equivalent conductivity. Using this criterion, $K_{NSF}$ is defined such that the net surface fluxes are identical in the heterogeneous and homogeneous media.

$$ \langle \mathbf{v} \rangle_{NSF} = \langle \mathbf{v} \rangle_{NSF} = \langle -K_{NSF} \nabla H \rangle_{NSF} \quad (75) $$

The brackets represent the averaging operator described in Equations (57)–(62) for the 2D case, and in Appendix A.5 for the 3D case. Recall that lower case **v** represents local microscale velocity (which can be highly variable in space), while upper case **V** represents macroscale velocity (which can be either constant or mildly variable in space). Caution is required with the $<...>$NSF averaging operator because it involves an integration of the normal fluxes on the boundary of the sample.

It is possible to express $K_{NSF}$ from Equation (75) after some manipulations. However, this leads to a complex averaging operator for the head gradient along the boundary of the domain. In fact, the problem lies not in this complex integration but in the fact that the spatial distribution of the head gradient along the boundary of the sample is itself a function of $K_{NSF}$. This means that in the general situation, one cannot manipulate Equation (75) to provide a closed form expression of $K_{NSF}$ without making an additional assumption.

When the hydraulic gradient $\mathbf{J} = -\nabla H$ is constant along the boundary of the domain of the homogenized porous block, Equation (75) has a closed form solution. This situation occurs, for example, when imposing linearly varying head, a constant vectorial flux, or periodic boundary conditions. When this is the case, the NSF formulas have been derived to ensure consistency, and we have $V_{NSF} = V$ under these conditions. We then obtain the following operational definition of the equivalent conductivity tensor from the NSF criterion with GRD, FLX, or PRD boundary conditions:

$$ \langle \mathbf{v} \rangle_{NSF} = K_{NSF} \mathbf{J} \quad (76) $$

When the head gradient is not constant in the homogeneous domain, one can use an estimation of the mean gradient within the domain and along the boundary to find an approximate solution to Equation (75). We propose, then, to use the volume average gradient criteria (VAG). This leads to the following operational definition for $K_{NSF}$.

$$ \langle \mathbf{v} \rangle_{NSF} \approx -K_{NSF} \langle \nabla h \rangle_{VAG} \quad (77) $$

Note that this last equation is identical to Equation (76) when using GRD, FLX, or PRD boundary conditions since the Volume Average Gradient in these cases is equal to the imposed gradient $-\mathbf{J}$. In the numerical experiments that will be described in Section 6, we will apply Equation (77) in all cases, that is, also in cases where the boundary conditions are not imposing a constant gradient in the homogeneous domain (PRM, for example), and we will illustrate numerically that the computation of $K_{NSF}$ using Equation (77) does not always ensure the equality Equation (75).

Like for the DIAG method, two numerical experiments are necessary to identify all the components of $K_{NSF}$ in 2D (or three experiments in 3D). Equation (77) is written for
each numerical experiment and assembled in a linear system where the unknown is the hydraulic conductivity tensor KNSF. This yields:

$$ \begin{pmatrix} K_{NSF}^{11} & K_{NSF}^{12} \\ K_{NSF}^{21} & K_{NSF}^{22} \end{pmatrix} \begin{pmatrix} < \partial_1 h >_{VAG}^1 & < \partial_1 h >_{VAG}^2 \\ < \partial_2 h >_{VAG}^1 & < \partial_2 h >_{VAG}^2 \end{pmatrix} = - \begin{pmatrix} < v_1 >_{NSF}^1 & < v_1 >_{NSF}^2 \\ < v_2 >_{NSF}^1 & < v_2 >_{NSF}^2 \end{pmatrix} \quad (78) $$

where ∂_i h stands for ∂h/∂x_i, and the superscript 1 or 2 associated with an average quantity is used to identify the flow experiment. The above equation can also be written in matrix notation:

$$ \mathbf{K}_{\text{NSF}} \mathbf{G}_{\text{VAG}} = -\mathbf{V}_{\text{NSF}} \quad (79) $$

In this equation, G<sub>VAG</sub> is a square matrix, whose columns contain the volume average gradient components computed for a given numerical experiment; V<sub>NSF</sub> is a square matrix as well, containing the corresponding Darcy velocities organized in a similar structure. The matrix representation is identical in 2D and 3D. If G<sub>VAG</sub> is invertible, the solution is obtained by multiplying both sides of the equation by G<sub>VAG</sub><sup>-1</sup>:

$$ \mathbf{K}_{\text{NSF}} = -\mathbf{V}_{\text{NSF}} \mathbf{G}_{\text{VAG}}^{-1} \quad (80) $$

In 2D, this explicit solution for the NSF equivalent permeability is simply:

$$ \mathbf{K}_{\text{NSF}} = \frac{1}{\delta_{\text{VAG}}} \begin{pmatrix} < v_1 >_{\text{NSF}}^1 & < v_1 >_{\text{NSF}}^2 \\ < v_2 >_{\text{NSF}}^1 & < v_2 >_{\text{NSF}}^2 \end{pmatrix} \begin{pmatrix} - < \partial_2 h >_{\text{VAG}}^2 & < \partial_1 h >_{\text{VAG}}^2 \\ < \partial_2 h >_{\text{VAG}}^1 & - < \partial_1 h >_{\text{VAG}}^2 \end{pmatrix} \quad (81) $$

with: δ<sub>VAG</sub> = < ∂<sub>1</sub>h ><sup>1</sup><sub>VAG</sub> < ∂<sub>2</sub>h ><sup>2</sup><sub>VAG</sub> - < ∂<sub>1</sub>h ><sup>2</sup><sub>VAG</sub> < ∂<sub>2</sub>h ><sup>1</sup><sub>VAG</sub>.

Note that Equation (81) contains all the components of the head gradient vectors for both numerical experiments. This allows computing the off-diagonal terms of the equivalent conductivity matrix K<sub>NSF</sub>. A similar procedure was employed by Kfoury et al. (2006).

To ensure that G<sub>VAG</sub><sup>-1</sup> is invertible, it is convenient to select the boundary conditions of the 2 in 2D or 3 in 3D numerical experiments such that the directions of the gradient are perpendicular. In particular, when constant gradient (GRD) or periodic conditions (PRD) are used, it is very convenient to use unit gradients aligned with the coordinates systems (<∇h><sup>i</sup><sub>VAG</sub> = -J<sup>i</sup> = -u<sub>i</sub> for i = 1, ..., 3), then G<sub>VAG</sub> becomes the identity matrix (G<sub>VAG</sub> = -I) and the previous equations become simply, under these unit gradient conditions:

$$ K_{\text{NSF}} = V_{\text{NSF}} = \begin{pmatrix} < v_1 >_{\text{NSF}}^1 & < v_1 >_{\text{NSF}}^2 & < v_1 >_{\text{NSF}}^3 \\ < v_2 >_{\text{NSF}}^1 & < v_2 >_{\text{NSF}}^2 & < v_2 >_{\text{NSF}}^3 \\ < v_3 >_{\text{NSF}}^1 & < v_3 >_{\text{NSF}}^2 & < v_3 >_{\text{NSF}}^3 \end{pmatrix} \quad (82) $$

As already explained earlier, the previous Equation (77) and all its simplifications just above ensure the equality of the NSF criterion between the homogeneous (homogenized) medium and the heterogeneous medium, when using either GRD or PRD boundary conditions (head gradient or permeametric conditions). However, Equation (77) remains an approximation in the general case.

## Equivalent K with the NSF Criterion for a 2D Fracture Network (Summary)

Briefly, in the case of a 2D fracture network, and for constant head gradient conditions (GRD), we have found that the equivalent permeability resulting from the previous NSF criterion can be expressed as:

$$ \mathbf{K}_{\text{NSF}} = \frac{1}{2L} \mathbf{P} \cdot (\mathbf{L}_{\text{FF}} - \mathbf{L}_{\text{IF}}^{\text{T}} \mathbf{L}_{\text{II}}^{-1} \mathbf{L}_{\text{IF}}) \cdot \mathbf{X} \quad (83) $$

In this expression (Equation (83)), K<sub>NSF</sub> represents the resulting 2 × 2 equivalent conductivity tensor (see comment on this further below). The other quantities are defined as follows:
■ Quantities such as $L_{II}$ designate node-node Laplacian matrices, weighted by the conductances $C$ of the fracture segments (links) relating the nodes (intersections). The discrete Laplacian matrix corresponds to the continuum Laplacian operator $div(K grad)$ but with discrete versions of the $div$ and $grad$ operators (they are expressed on the graph representing the fracture network). The conductivities $K$ are replaced by conductances $C$ based on either Darcy or Poiseuille's law for flow through each link. The Laplacian matrix to be inverted in Equation (83) is the internal node-node matrix $L_{II}$ of size $N_I \times N_I$. Other auxiliary node-node Laplacian matrices appear also; they are indexed depending on whether they concern internal nodes (I) and/or boundary nodes (F): thus, $L_{II}$ is the Laplacian matrix of size $N_I \times N_I$ relating all internal nodes, while $L_{IF}$ of size $N_I \times N_F$ is the Laplacian matrix relating internal nodes to boundary nodes, and $L_{FF}$ of size $N_F \times N_F$ is the Laplacian relating all boundary nodes between them (if there does not exist a single link between any pair of boundary nodes, then $L_{FF} = 0$).

■ $P$ represents a geometric matrix of size $2 \times N_F$ containing the two components of the unit vector normal to the block at the $N_F$ boundary nodes (which are submitted to linearly varying head for the present case of gradient condition “GRD”);

■ $X$ represents a coordinate matrix of size $N_F \times 2$ with $N_F$ lines (number of boundary nodes) and two-columns $[x \ y]$ containing the coordinates of the $N_F$ boundary nodes;

■ $L$ is a scalar representing the size of the 2D fractured domain (here a square block); note: the domain size “L” should not be confused with the Laplacian matrix “L”.

With these definitions, the reader can check that the resulting equivalent conductivity $K_{NSF}$ is indeed an indexed $2 \times 2$ matrix $(K_{ij})_{NSF}$; furthermore, analyzing the algebraic nature of the expression of Equation (83), it can be shown that the resulting $(K_{ij})_{NSF}$ behaves like a second rank tensor. Properties of this equivalent tensor, and those obtained by applying other criteria, have been analyzed by us [61,62]) and will be presented in more detail elsewhere in future (e.g., [97]).

For more details on the algebra of graphs, see Gil Strang's seminal paper on equilibrium problems in hydraulics and structural mechanics [98], where graph structures and algebraic approaches are emphasized. His ideas were also presented more recently in a book ([99], Chap. 2: A framework for applied mathematics therein).

### 4.2.3. K from Vectorial Surface Flux (VSF)

In a similar manner, we can define the equivalent conductivity tensor $K_{VSF}$ such that the Vectorial Surface Flux is identical on both the homogeneous and heterogeneous media. Here, the averaging operator is simpler than for the VSF expression and can be expressed as surface integrals on the local and homogenized Darcy velocities (i.e., the microscale and macroscale velocities).

$$ \langle \mathbf{v} \rangle_{\text{VSF}} = \langle \mathbf{v} \rangle_{\text{VSF}} \quad \text{or} \quad \frac{1}{\Gamma} \int_{\Gamma} \mathbf{v} \, d\gamma = \frac{1}{\Gamma} \int_{\Gamma} \mathbf{V} \, d\gamma \qquad (84) $$

Replacing the large-scale specific discharge $V$ by the upscaled Darcy's law with a constant tensor $K_{VSF}$ leads to the following definition of the equivalent conductivity for continuous or composite media:

$$ \langle \mathbf{v} \rangle_{\text{VSF}} = -K_{\text{VSF}} \int_{\Gamma} \nabla H \, d\gamma \qquad (85) $$

Like for the NSF method, the average value of the gradient on the frontier of the block of the homogeneous medium is not known in general as it is a function of $K_{VSF}$, the geometry of the block, and the boundary conditions. If the boundary conditions are such that the gradient is constant in the homogeneous medium, then this average quantity is known and independent from the value of $K_{VSF}$. It is equal to the imposed gradient. If this
is not the case, we use the volume average gradient on the heterogeneous medium as a reasonable approximation.

$$ \langle \mathbf{v} \rangle_{\text{VSF}} \approx -K_{\text{VSF}} \langle \nabla h \rangle_{\text{VAG}} \quad (86) $$

To compute $K_{VSF}$, Equation (86) is used twice in 2D or three times in 3D with different numerical experiments in order to identify all the components of the conductivity tensor with an expression analogous to Equation (80).

$$ \mathbf{K}_{\text{VSF}} = -\mathbf{V}_{\text{VSF}} \cdot \mathbf{G}_{\text{VAG}}^{-1} \quad (87) $$

The only difference with Equation (80) is that the averaging technique used to compute the mean velocity is here the Vectorial Surface Flux while it was the Net Surface Flux in the previous method. Again, if linearly varying head or periodic boundary conditions with unit gradients aligned with the system of coordinates, one can use a similar expression as Equation (82) with the VSF averages.

### Equivalent K with the VSF Criterion for a 2D Fracture Network (Summary)

Briefly, in the case of a 2D fracture network, and for constant head gradient conditions (GRD), we have found that the equivalent permeability resulting from the previous NSF criterion can be expressed as:

$$ K_{\text{VSF}} = \frac{1}{2L} \mathbf{P} \cdot (\mathbf{L}_{\text{FF}} - \mathbf{L}_{\text{IF}}^{\text{T}} \mathbf{L}_{\text{II}}^{-1} \mathbf{L}_{\text{IF}}) \cdot \mathbf{X} \quad (88) $$

where (briefly): $U$ is a $2 \times N_D$ matrix indicating the direction vectors of the $N_D$ "boundary" fracture links that intersect the boundary of the domain; $S_{DF}$ and $S_{DI}$ are submatrices that relate (through conductances) those $N_D$ links to (respectively) the $N_F$ boundary nodes (F) and the $N_I$ interior nodes (I); $L_{II}$ is the square Laplacian matrix of size $N_I \times N_I$ over all interior nodes; $L_{IF}$ is the $N_I \times N_F$ Laplacian matrix relating interior to boundary nodes; and $X$ is an $N_F \times 2$ coordinate matrix containing the two coordinates of the $N_F$ boundary nodes. Note that matrices $S$ are of the form $S = C A^T$, where $C$ is the diagonal conductance matrix and $A$ is the node-node adjacency matrix. The resulting equivalent permeability $K_{VSF}$ is given explicitly by Equation (88), in an algebraic fashion. The main computational task is (as before) the inversion of the square Laplacian matrix $L_{II}$.

#### 4.2.4. Volume Average Flux (VAF)

A third possibility is to define the equivalent conductivity tensor $K_{VAF}$ such that the Volume Average Flux (VAF) is identical on both the homogeneous and heterogeneous media. Here, the averaging operator is a volume integral, although we have already seen that it can be replaced by a surface integral without any loss of generality. The volume average flux criterion is:

$$ \frac{1}{\Omega} \int_{\Omega} v d\omega = \frac{1}{\Omega} \int_{\Omega} V d\omega \quad \text{or} \quad \langle \mathbf{v} \rangle_{\text{VAF}} = \langle \mathbf{v} \rangle_{\text{VAF}} \quad (89) $$

Replacing the macroscale Darcy velocity $V$ by the upscaled Darcy's law leads to the following definition of the upscaled VAF conductivity $K_{VAF}$:

$$ \langle \mathbf{v} \rangle_{\text{VAF}} = -K_{\text{VAF}} \int_{\Omega} \nabla H d\omega \quad (90) $$

As before, $K_{VAF}$ is implicitly in the volume integral because $H$ is a function of the geometry of the block, the boundary conditions, and $K_{VAF}$. When constant gradients (GRD) or periodic boundary (PRD) conditions are used, then the Darcy velocity and the head gradient are constant in the domain, whatever is its geometry, and the volume average gradient is directly known. In the other cases, we will invoke an additional criterion, which

Erratum for Eq.(88) and the text above it; the correct text and equation are given below:

Equivalent K with the VSF criterion for a 2D fracture network (summary)

In the case of a 2D fracture network, and for constant head gradient conditions (GRD), the equivalent permeability resulting from the VSF criterion is of the form:

$$ K_{\text{VSF}} = \frac{1}{2L} \mathbf{U} \cdot (S_{DF} - S_{DI}A^{-1}A_{IF}) \cdot \mathbf{X} \quad \text{Eq. 88} $$
is that the volume averaged head gradients in both the heterogeneous and homogeneous domains should be equal:

$$ \frac{1}{\Omega} \int_{\Omega} \nabla h \, d\omega = \frac{1}{\Omega} \int_{\Omega} \nabla H d\omega \qquad (91) $$

Using both the VAF and VAG criteria leads to the operational definition of the equivalent conductivity used by many authors since Rubin and Gomez-Hernandez (1990):

$$ \langle \mathbf{v} \rangle_{\text{VAF}} = -\mathbf{K}_{\text{VAF}} \int_{\Omega} \nabla h \, d\omega \qquad (92) $$

Here, it is important to note that the reciprocal may not be true, computing $\mathbf{K}_{\text{VAF}}$ from Equation (92) does not necessarily imply that the VAF and VAG criteria are both verified.

As for $\mathbf{K}_{\text{VAF}}$ or $\mathbf{K}_{\text{NSF}}$, $\mathbf{K}_{\text{VAF}}$ is computed numerically by solving two numerical problems in 2D or three in 3D and by applying Equation (92) for each numerical experiment. All those equations are assembled to identify all the components of the conductivity tensor.

$$ \mathbf{K}_{\text{VAF}} = -\mathbf{V}_{\text{VAF}} \cdot \mathbf{G}_{\text{VAG}}^{-1} \qquad (93) $$

with $\mathbf{G}_{\text{VAG}}$ and $\mathbf{V}_{\text{VAF}}$ being the matrices containing the compilation of the volume average fluxes and gradients computed from the results of two numerical experiments with perpendicular gradients as explained in detail earlier (see Equations (78) and (81)).

As before, this equation can be simplified for linearly varying head or periodic conditions with unit gradients.

$$ K_{\text{VAF}} = V_{\text{VAF}} = \begin{pmatrix} < v_1 >_{\text{VAF}}^1 & < v_1 >_{\text{VAF}}^2 & < v_1 >_{\text{VAF}}^3 \\ < v_2 >_{\text{VAF}}^1 & < v_2 >_{\text{VAF}}^2 & < v_2 >_{\text{VAF}}^3 \\ < v_3 >_{\text{VAF}}^1 & < v_3 >_{\text{VAF}}^2 & < v_3 >_{\text{VAF}}^3 \end{pmatrix} \qquad (94) $$

Equivalent K with the VAF Criterion for a 2D Fracture Network (Summary)

Briefly, in the case of a 2D fracture network, and for constant head gradient conditions (GRD), the equivalent permeability resulting from the VAF criterion is of the form:

$$ \mathbf{K}_{\text{VAF}} = \frac{1}{L^2} \mathbf{X}^{\text{T}} \cdot (\mathbf{L}_{\text{FF}} - \mathbf{L}_{\text{IF}}^{\text{T}} \mathbf{L}_{\text{II}}^{-1} \mathbf{L}_{\text{IF}}) \cdot \mathbf{X} \qquad (95) $$

where the various matrices are defined as in Equation (83), with Laplacian matrices $\mathbf{L}_{\text{FF}}$ (boundary nodes to boundary nodes), $\mathbf{L}_{\text{IF}}$ (interior nodes to boundary nodes), $\mathbf{L}_{\text{II}}$ (interior nodes to interior nodes), and $\mathbf{X}$ is the $N_F \times 2$ rectangular matrix of coordinate positions of the $N_F$ boundary nodes. Again, the 2D block is assumed here to be a square of side "L" (not to be confused with a Laplacian matrix).

### 4.2.5. Equivalent K with the Volume Average Gradient (VAG) Criterion

We have seen that in the previous definitions of the equivalent conductivity, we often had to use the volume average gradient as complementary information. In addition, when constant flux boundary conditions are used, all the previous criteria are directly imposed by the boundary conditions and cannot be used to define $\mathbf{K}$. In this case, one has to use the Volume Average Gradient to define and compute the equivalent conductivity $\mathbf{K}_{\text{VAG}}$.

$$ \frac{1}{\Omega} \int_{\Omega} \nabla h d\omega = \frac{1}{\Omega} \int_{\Omega} \nabla H d\omega \qquad (96) $$

Using the large-scale Darcy's law, and taking $\mathbf{K}_{\text{VAG}}$ out of the averaging operator, one obtains:

$$ K_{\text{VAG}} \int_{\Omega} \nabla h \, d\omega = - \int_{\Omega} V d\omega \qquad (97) $$
Like for the previous cases, one can fully identify $K_{VAG}$ by running three different flow experiments. When using constant flux (FLX) boundary conditions, the right-hand side of Equation (97) is imposed by the boundary conditions. If we assume that the two numerical experiments are designed with unit vectorial fluxes oriented along the $x_1$ and $x_2$ directions, we get the following expression for $K_{VAG}$.

$$ K_{VAG} = -G_{VAG}^{-1} = - \begin{pmatrix} < \partial_1 h >_{VAG}^1 & < \partial_1 h >_{VAG}^2 & < \partial_1 h >_{VAG}^3 \\ < \partial_2 h >_{VAG}^1 & < \partial_2 h >_{VAG}^2 & < \partial_2 h >_{VAG}^3 \\ < \partial_3 h >_{VAG}^1 & < \partial_3 h >_{VAG}^2 & < \partial_3 h >_{VAG}^3 \end{pmatrix}^{-1} \quad (98) $$

### 4.2.6. Equivalent K from the Total Dissipated Power (TDP) Criterion

Some authors such as [30,94,96,100,101] used the Total Dissipated Power to define the equivalent hydraulic conductivity. As in the previous definitions, the principle here is to state that the total dissipated power must be identical in the heterogeneous medium and in the corresponding homogeneous (homogenized) medium:

$$ - \rho g \int_{\Omega} (\mathbf{v} \cdot \nabla h) \, d\omega = - \rho g \int_{\Omega} (\mathbf{V} \cdot \nabla H) \, d\omega \quad (99) $$

If we replace, in the right-hand side, the block scale or macroscale $V$ by the macroscale Darcy phenomenological law, we then obtain the following definition of the equivalent conductivity matrix $K_{TDP}$:

$$ \int_{\Omega} (\mathbf{v} \cdot \nabla h) \, d\omega = \int_{\Omega} (K_{TDP} \nabla H \cdot \nabla H) \, d\omega \quad (100) $$

This criterion has been used very successfully for deriving bounds for the equivalent conductivity [43,94]. However, because this criterion leads to a scalar equation, it is a priori less convenient than a vectorial criterion if the objective is to calculate a non-diagonal equivalent conductivity matrix (or more generally, the equivalent conductivity tensor). In particular, because each numerical experiment will provide only one equation, this criterion (Equation (100)) will require running more flow simulations than the other techniques.

To avoid that problem, Bøe [96] proposes to split the implicit sums resulting from the product of the head gradient with all the components of the Darcy velocity in the criterion of Equation (100). This results in a set of nine sufficient criteria involving a basis of three numerical flow problems in perpendicular directions whose solutions are $h^1$ to $h^3$. Assuming that the boundary conditions are periodic, Bøe [96] proposes an operational definition, which can be written with our own notations as follows.

$$ K_{ir} = \frac{< v_i^r >_{VAF}}{J_i} \quad (101) $$

This expression (Equation (101)) is a sufficient condition to satisfy the TDP criterion. It was also suggested by Farmer [35] as "a natural" definition of the equivalent conductivity without stating a specific equivalence criterion. Both Bøe [96], assuming a diagonal local conductivity tensor and periodic boundary conditions, and Farmer [35], for a symmetric but not necessarily diagonal local conductivity tensor and linearly varying head, show that the previous equation can be written in a more complex form.

$$ K_{rs} = \frac{1}{J_r J_s} \int_{\Omega} \left( \sum_i \sum_j k_{ij} \frac{\partial h^r}{x_j} \frac{\partial h^s}{x_i} \right) d\omega \quad (102) $$

Equation (102) is not used to compute $K_{TDP}$, but it directly shows the symmetry of $K_{TDP}$ obtained from Equation (101) if one uses linearly varying head or periodic boundary conditions. If the intensity of the gradient is $J_i = 1$ in each numerical experiment, Equation (101) is identical to Equation (94). Furthermore, the resulting tensor is strictly
the same as the one computed by the VAF technique when using GRD or PRD boundary conditions. Therefore, in the remainder of this paper, we will not consider that approach as a distinct approach to compute **K**. Instead, because the TDP is an important criterion, we will test whether all the other definitions of the equivalent conductivity tensor will satisfy the total dissipated power.

### 4.3. Dual Permeability/Resistivity Matching

We end up with a brief presentation of the dual permeability/resistivity matching procedure for determining both the equivalent permeability and equivalent resistivity. This procedure will not be implemented numerically to the various block samples studied in Section 6 (some of them random, others deterministic), but it is worth mentioning, as it has proven useful in the context of upscaling flow in randomly heterogeneous media via analytical perturbation methods (e.g., [75,102,103]).

This approach is based on the local head governing PDE $\nabla \cdot (\mathbf{k}(\mathbf{x}) \nabla h(\mathbf{x})) = 0$, but uses it in two different ways to obtain an upscaled Darcy law (two “dual” approaches):

* upscaling is performed in terms of the local hydraulic conductivity **k** (based on Darcy's law $\mathbf{v} = -\mathbf{k}(\mathbf{x}) \nabla h$), leading to an equivalent macroscale conductivity **K**, and ...

* upscaling is also conducted in terms of hydraulic resistivity $\mathbf{r} = \mathbf{k}^{-1}$ (based on the reciprocal formulation of Darcy's law $-\nabla h = \mathbf{r}(\mathbf{x}) \cdot \mathbf{v}$), leading to an equivalent macroscale resistivity **R**.

The term “conjugate” or “dual” was introduced in Ababou ([1], Chap.4) to describe the particular invariance properties of 2D heterogeneous continua with equilibrated $\ln k(\mathbf{x})$ distributions, under the transformation $k \to k^{-1}$, exploiting and generalizing previous results by Shvidler [73] and Matheron [30]. The term “dual upscaling” was later applied to the tensorial macro-resistivity $R_{ij}$ obtained by upscaling $-\nabla h = \mathbf{r}(\mathbf{x}) \cdot \mathbf{v}$ in any number of spatial dimensions, where $\mathbf{r}(\mathbf{x})$ is the local resistivity (e.g., [102] for single phase flow in a 3D random medium). This “dual” approach leads to a comparison of Kij and Rij based on algebraic positivity and reciprocity upscaling criteria: $[R_{ij}] > 0$; $[K_{ij}] > 0$; $[R_{ij}] = [K_{ij}]^{-1}$.

These criteria were presented in Fadili and Ababou [75] for a random continuum with spatially correlated log-parameters, first for single phase flow, and secondly for immiscible two-phase flow based on Darcy-Muskat equations, with nonlinear retention curves and two-phase permeability curves. This permeability/resistivity approach had also been presented earlier, for single-phase Darcy flow in randomly heterogeneous media (e.g., [102,103]). In summary, the main result for the case of a statistical Darcian continuum with random field Gaussian log-conductivity $F = \ln K(x, y, z)$ is of the form:

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

$$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}$$

See also around Equation (101) in Appendix A.8. This formulation 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 (for more details, see Appendix A.8). The permeability/resistivity matching approach could be considered in the future for application to finite domain upscaling, that is, matching macroscale block-permeability and block-resistivity.

## 5. Explicit Expressions of Equivalent Conductivity Tensors for Layered Media

In this section, we develop expressions of equivalent conductivity tensors for a few special types of heterogeneous media or fracture networks with prescribed boundary conditions.
This section shows that piecewise linear head conditions lead (analytically) to a tensorial Kij. Two methods were developed, one by Renard ([80], Annex B therein), and the other was presented by Cañamón et al. [82,83] and Ababou et al. [84]. These methods yield the same types of results for slightly different layered configurations. In the latter method, the microscale velocity and head gradient fields are piecewise constant, and the microscale head field is piecewise linear, both on the boundary and inside the layered domain. Furthermore, the microscale and macroscale head fields (h and H) do not obey the same boundary conditions, unlike the gradient condition “GRD” with linearly varying head condition, presented in Section 2.4. The present approach with piecewise linear head conditions can be described as follows:

* Conditions: immersion in a piecewise constant gradient (GRD') or flux (FLX')

* Criteria: Volume Average Gradient and Volume Average Flux (VAG+VAF)

## 5.1. Exact Flow Solution and Equivalent Permeability for a Three-Layered Fractured Medium under Piecewise Linear Head Conditions (GRD')

We present here a special analysis that produces an exact tensorial Kij for multi-layered 2D or 3D prismatic domains submitted to specific piecewise linear head distributions on the boundary (GRD'). Among the various upscaling methods analyzed in this paper, these particular head boundary conditions were not fully investigated in the previous section; however, this case is still relevant because the simple flow solution obtained with these conditions on a layered block leads to an exact evaluation of the equivalent block permeability tensor Kij, which can be used as a reference for comparisons with other methods. Cañamón [82], Cañamón et al. [83,104]), Ababou et al. [61], applied the method to a composite matrix/fracture “unit block”, as the first step in a sequential upscaling procedure for 3D fractured porous media with a permeable isotropic matrix. The fractures can behave either isotropically (assuming the fractures are similar to thin layers of coarse porous material) or anisotropically (assuming Poiseuille flow parallel to fracture planes). The upscaling presented below is also applicable to N layers instead of a fractured porous matrix.

Here, let us focus on the case of a unit block comprising three “layers”: an inner “fracture” or “coarse” layer (with isotropic Darcy permeability k<sub>F</sub> or with directional Poiseuille permeability k<sub>F</sub>''), bordered by 2 “matrix” layers, each with isotropic Darcy permeability k<sub>M</sub> < k<sub>F</sub>. The geometry and boundary conditions for the “unit fractured block” problem are illustrated in Figure 11 for a rectangular box-shaped domain, although more general prismatic shapes yield similar results.

To upscale the flow on this geometry, we first look for a local velocity distribution

$$ (v(x))_{\Omega_k} = (v_x, v_y, v_z)_{\Omega_k} $$

such that, on each sub-domain Ω<sub>k</sub> (k = A,B,C), except at the material interfaces Γ<sub>I(k)</sub>, we have:

$$ \left. \begin{array}{l} \mathrm{div}(v) = \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} + \frac{\partial v_z}{\partial z} = 0 \\ v = k \quad j = -k \operatorname{grad}(h) \end{array} \right\} \forall \mathbf{x} \in \Omega_k \ (x \notin \Gamma_{I(k)}) \qquad (104) $$

submitted to piecewise linear head boundary conditions of as follows:

$$ h(x, y, z) = -j_0^{(k)} \bullet \mathbf{x} = -(j_{0x}^{(k)} x + j_{0y}^{(k)} y + j_{0z}^{(k)} z) \quad \forall (x, y, z) \in \Gamma_F^{(k)} \qquad (105) $$

**Notations.** All lower-case quantities are local microscale quantities. They are generally non constant in space. For instance, the hydraulic gradient j<sub>0</sub> is imposed (subscript 0) but not necessarily constant in space (lower case j). In the present case, because j<sub>0</sub> is piecewise constant, we use the notation j<sub>0</sub><sup>(k)</sup>, where index (k) describes the sub-domains.
<figure><img src="image_14.png" /><figcaption>**Figure 11. (a-left):** Schematic of a "unit" block Ω containing a single "fracture" Ω<sub>F</sub> sandwiched between two layers of porous matrix Ω<sub>M</sub>. **(a-right):** The fractured matrix block, with explicit indication of fracture aperture and matrix layers thicknesses. **(b):** Two-dimensional view of the piecewise linear head conditions, corresponding to piecewise constant gradient **J** along the boundaries (GRD').</figcaption></figure>

The above condition is applied separately on each piece of the boundary (on each subdomain k = A,B,C). Note that the boundary conditions are chosen here such that h(x,y,z) is of the lowest possible order: h(x, y, z) C<sup>0</sup>(R<sup>3</sup>) (space of continuous functions). The head field is not necessarily differentiable, or possibly, only differentiable along some directions. For this reason, the final upscaling result, although exact for these boundary conditions, is a "low order" equivalent conductivity tensor (from that point of view).

The next step is to implement mass flux continuity and pressure continuity along all "layer/layer" interfaces Γ<sub>I(k)</sub>, that is ∀(x, y, z) ∈ Γ<sub>I(k)</sub>, including at boundary points where interfaces meet boundaries, that is ∀x ∈ Γ<sub>I(k)</sub>. For instance, in the case of a horizontal interface, the continuity relations are.

$$
\left\{
\begin{array}{l}
(\mathbf{v} \cdot \mathbf{n})^{+} = (\mathbf{v} \cdot \mathbf{n})^{-} \Rightarrow \mathbf{v}_{z}^{+} = \mathbf{v}_{z}^{-} \\
h^{+}(\mathbf{x}) = h^{-}(\mathbf{x}) \Rightarrow \left(\frac{\partial h}{\partial x}\right)^{+} = \left(\frac{\partial h}{\partial x}\right)^{-} \Leftrightarrow \mathbf{j}_{x}^{+} = \mathbf{j}_{x}^{-}
\end{array}
\right.
\qquad (106)
$$

Notations for interfaces: Γ<sub>I(k)</sub> designates all interfaces pertaining to sub-domain Ω<sub>k</sub>; on the other hand, a notation such as Γ<sub>A/B</sub> can be used to denote the interface formed by the intersection of the two sub-domains Ω<sub>A</sub> and Ω<sub>B</sub> (see Figure 11).

After some manipulations, which are not detailed here, and expressing the piecewise
constant gradient j₀⁽ᵏ⁾ in terms of a constant reference gradient denoted J₀, a unique set of
piecewise linear head conditions (and gradients) is obtained that satisfies the continuity
conditions, not only inside the domain but also on its boundary. The reader can verify this by direct insertion of the following equation:

$$ \mathbf{j}(\mathbf{x}) = \begin{vmatrix} \mathbf{j}_M & \text{if } \mathbf{x} \in \Omega_M \\ \mathbf{j}_F & \text{if } \mathbf{x} \in \Omega_F \end{vmatrix} \text{ with } \mathbf{j}_M = \mathbf{J}_0; \mathbf{j}_F = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & K_M/K_F \end{pmatrix} \bullet \mathbf{J}_0 \quad (107) $$

These boundary conditions finally produce a piecewise constant velocity vector $\mathbf{v}(\mathbf{x})$, which can be controlled by controlling $\mathbf{J}_0$, and which satisfies mass balance and continuity conditions throughout the composite domain ($\Omega$):

$$ \mathbf{v}(\mathbf{x}) = \begin{vmatrix} \mathbf{v}_M & \text{if } \mathbf{x} \in \Omega_M \\ \mathbf{v}_F & \text{if } \mathbf{x} \in \Omega_F \end{vmatrix} \quad (108) $$

$$ \text{where } \mathbf{v}_M = \begin{pmatrix} k_M & 0 & 0 \\ 0 & k_M & 0 \\ 0 & 0 & k_M \end{pmatrix} \bullet \mathbf{J}_0; \quad \mathbf{v}_F = \begin{pmatrix} k_F^{\parallel} & 0 & 0 \\ 0 & k_F^{\parallel} & 0 \\ 0 & 0 & k_M \end{pmatrix} \bullet \mathbf{J}_0 \quad (109) $$

Using this flow field, the upscaled conductivity for this special problem is obtained below.

In order to upscale this local flow problem, the average flux and gradient vectors are computed at block scale with one of the available methods, given the exact “local” flow field given by Equations (108) and (109). Thus, using the VAF method for the block scale flux and the VAG method for the block scale gradient, we obtain respectively (There are in fact several ways to proceed with this average. Here we consider the Volume Averaged Flux (VAF). However, in Appendix IX of Canamon [82], the approach based on Vectorial Surface Flux (VSF) is also presented (which may be useful for comparison purposes).):

$$ \mathbf{V} = \frac{\int_{\Omega} \mathbf{v} \, d\omega}{\int_{\Omega} d\omega} \Rightarrow \mathbf{V} = \frac{(V_M)_{\Omega}}{V_{\Omega}} \cdot \mathbf{v}_M + \frac{(V_F)_{\Omega}}{V_{\Omega}} \cdot \mathbf{v}_F \quad (110) $$

$$ \mathbf{J} = \frac{\int_{\Omega} \mathbf{j} \, d\omega}{\int_{\Omega} d\omega} \Rightarrow \mathbf{J} = \frac{(V_M)_{\Omega}}{V_{\Omega}} \cdot \mathbf{j}_M + \frac{(V_F)_{\Omega}}{V_{\Omega}} \cdot \mathbf{j}_F \quad (111) $$

Denoting $\phi$ the volumetric fraction of the inner “fracture” layer within the block, and $(1 - \phi)$ the volumetric fraction of the matrix, we obtain the block scale (VAF) flux vector as:

$$ \mathbf{V} = \begin{pmatrix} (1 - \phi) \cdot k_M + \phi \cdot k_F^{\parallel} & 0 & 0 \\ 0 & (1 - \phi) \cdot k_M + \phi \cdot k_F^{\parallel} & 0 \\ 0 & 0 & k_M \end{pmatrix} \bullet \mathbf{J}_0 \quad (112) $$

and inserting the block scale (VAG) gradient $\mathbf{J}$, we obtain a “Darcian” linear tensorial relation between the macroscopic flux and gradient vectors, as follows:

$$ \mathbf{V} = \begin{pmatrix} (1 - \phi) \cdot k_M + \phi \cdot k_F^{\parallel} & 0 & 0 \\ 0 & (1 - \phi) \cdot k_M + \phi \cdot k_F^{\parallel} & 0 \\ 0 & 0 & \frac{1}{\frac{(1-\phi)}{k_M} + \frac{\phi}{k_F^{\perp}}} \end{pmatrix} \bullet \mathbf{J}_0 = \mathbf{K} \mathbf{J} \quad (113) $$

Finally, this yields the equivalent tensorial conductivity $\mathbf{K}$ of the block in the block’s reference frame:

$$ \mathbf{K} = \begin{pmatrix} K_{11} & 0 & 0 \\ 0 & K_{22} & 0 \\ 0 & 0 & K_{33} \end{pmatrix} = K_{ij} \cdot \delta_{ij}; \quad \left\{ \begin{array}{l} K_{11} = (1 - \phi) \cdot k_M + \phi \cdot k_F^{\parallel} \\ K_{22} = (1 - \phi) \cdot k_M + \phi \cdot k_F^{\parallel} \\ K_{33} = \frac{1}{\frac{(1-\phi)}{k_M} + \frac{\phi}{k_F^{\perp}}} \end{array} \right. \quad (114) $$
where the resulting principal conductivities (Kii*) can be interpreted as weighted arithmetic and harmonic means. More generally, if the 3D block is oriented at angle with respect to the reference frame, we obtain:

$$ (K_{ij})_{\Omega} = (\delta_{ij} - n_i n_j) \cdot K_A + n_i n_j \cdot K_H \quad (115) $$

where: $(n_i, i = 1,2,3)$ are the components of the unit vector **n** normal to the layers; $K_A$ and $K_H$ are the weighted arithmetic and harmonic mean conductivities (defined above).

Specialization: Equivalent Kij for a Single Fracture Block with Impervious Matrix

Finally, it is interesting to specialize the above result to the case of a single fracture traversing an impervious block with zero matrix permeability; this is a very basic building piece of any “fracture network”. Thus, we set $k_M = 0$, and we obtain as a result a tensorial block conductivity of the form:

$$ (K_{ij})_{\Omega} = \phi \cdot k_F^{\parallel} \cdot (\delta_{ij} - n_i n_j) \quad (116) $$

where $\phi$ is simply the volumetric fraction of the single “fracture” within the impervious block. For a Poiseuille fracture, the conductivity $k_F^{\parallel}$ in Equation (116) is the Poiseuille conductivity parallel to the fracture plane (law of cubic aperture). However, if the fracture is made up of (or filled with) a coarse isotropic porous material, then we set $k_F^{\parallel}$ in Equation (116) to represent the local isotropic Darcy conductivity inside the coarse “fracture” medium.

## 5.2. Exact Flow Solution for Two Anisotropic Layers under Piecewise Linear Head Conditions (GRD')

Renard ([80], Appendix B) developed a similar 3D solution under piecewise linear head conditions for a two-layer Darcian medium with fully non-diagonal permeability tensors, one in each layer. Thus, let us consider a composite medium made of two adjacent homogeneous 3D blocks, as shown in Figure 12.

The two blocks have the shape of a rectangular cuboid, and they share a common face. They both have a constant hydraulic conductivity tensor $k_1$ and $k_2$. The aim of this section is to derive an analytical expression for the equivalent conductivity tensor K that represents globally the composite block under the assumptions (1) that the local conductivity tensors $k_1$ and $k_2$ are symmetrical definite positive and constant in each block and (2) that the local scale-specific discharge rates $v_1(x)$ and $v_2(x)$ are different but spatially constant in each block (the streamlines are straight lines in each block). Under these assumptions, the head gradients are constant (but different) within each domain, and one can derive an exact analytical solution. A gradient in three dimensions can be defined by the values of four heads at four locations. To define these gradients and the continuity of head values at the boundary between the blocks, we introduce head values at some corners of the domain (see Figure 12). The values that are not named in the figure depend uniquely on the ones that are defined, and it is therefore not needed to introduce them.
<figure><img src="image_15.png" /><figcaption>Figure 12. Composite medium made of two rectangular cuboid blocks. Definition of the geometry and of the boundary conditions.</figcaption></figure>

To obtain the solution, the first step consists in defining the boundary conditions, which ensures that assumption (2) is honoured. To get the most general solution under this framework, we assume in addition that none of the head gradient component is equal to zero. With the geometry and notations introduced in Figure 12, the expression for the head gradient $(\nabla h)_1$ in block 1 (left side) is:

$$ (\nabla h)_1 = \begin{pmatrix} \frac{h_1-h}{a} \\ \frac{h_0-h_1}{h_4-h_1} \\ \frac{c}{d} \end{pmatrix} \qquad (117) $$

and for the head gradient $(\nabla h)_1$ in block 2 (right side), we have:

$$ (\nabla h)_2 = \begin{pmatrix} \frac{h-h_2}{b} \\ \frac{h_3-h_2}{h_5-h_3} \\ \frac{c}{d} \end{pmatrix} \qquad (118) $$

The continuity of the heads at the boundary between the two blocks implies the equality of the components of the gradient in the y and z (2 and 3) directions for the two blocks, and therefore we observe that we have:

$$ h_3 = h_2 + h_0 - h_1 \qquad (119) $$

$$ h_5 = h_3 + h_4 - h_1 \qquad (120) $$

To identify the head value noted h in Figure 12, we write the continuity of the local scale flux normal to the interface between the two blocks. This implies that the x component of the local scale Darcy velocity in the two blocks must be equal.

$$ v_1^1 = v_1^2 \qquad (121) $$

In the notation $v_i^j$, the subscript i represents the ith component of the velocity vector, and the superscript j allows to identify the block; $v_1^2$ is the x component of the local Darcy
velocity in block 2. By applying the local scale flow equation on both blocks, one obtains the equality:

$$k_1^{11} \left( \frac{h_1 - h}{a} \right) + k_1^{12} \left( \frac{h_0 - h_1}{c} \right) + k_1^{13} \left( \frac{h_4 - h_1}{d} \right) = k_2^{11} \left( \frac{h - h_2}{b} \right) + k_2^{12} \left( \frac{h_0 - h_1}{c} \right) + k_2^{13} \left( \frac{h_4 - h_1}{a} \right) \quad (122)$$

from which the unknown h can be isolated. Hence:

$$h = \frac{b k_1^{11} h_1 + a k_2^{11} h_2 + ab \left[ (k_1^{12} - k_2^{12}) \left( \frac{h_0 - h_1}{c} \right) + (k_1^{13} - k_2^{13}) \left( \frac{h_4 - h_1}{d} \right) \right]}{b k_1^{11} + a k_2^{11}} \quad (123)$$

Then, we define the equivalent conductivity K as the conductivity of a homogeneous block such that

$$\begin{pmatrix} V_1 \\ V_2 \\ V_3 \end{pmatrix} = - \begin{pmatrix} K^{11} & K^{12} & K^{13} \\ K^{21} & K^{22} & K^{23} \\ K^{31} & K^{32} & K^{33} \end{pmatrix} \begin{pmatrix} \frac{h_1 - h_2}{a+b} \\ \frac{h_0 - h_1}{h_4 - h_1} \\ \frac{c}{d} \end{pmatrix} \quad (124)$$

and such that the upscaled fluxes that traverse all the sides of the homogeneous medium should be identical to the ones traversing the composite heterogeneous medium. Let us first express the component of the Darcy velocity $V_1$ in the first direction (or rather, its opposite $-V_1$). It is identical in both blocks and equal to:

$$-V_1 = k_1^{11} \left( \frac{h_1 - h}{a} \right) + k_1^{12} \left( \frac{h_0 - h_1}{c} \right) + k_1^{13} \left( \frac{h_4 - h_1}{d} \right) \quad (125)$$

If we replace h by the result of our previous computation (Equation (123)), we obtain the expression

$$-V_1 = \left[ \frac{(a+b) k_1^{11} k_2^{11}}{b k_1^{11} + a k_2^{11}} \right] \left( \frac{h_1 - h_2}{a} \right) + \left[ \frac{b k_1^{11} k_2^{12} + a k_1^{11} k_2^{12}}{b k_1^{11} + a k_2^{11}} \right] \left( \frac{h_0 - h_1}{c} \right) \\ + \left[ \frac{b k_1^{11} k_2^{13} + a k_2^{11} k_1^{13}}{b k_1^{11} + a k_2^{11}} \right] \left( \frac{h_4 - h_1}{d} \right) \quad (126)$$

In the equivalent homogeneous medium, the same flux is expressed as

$$-V_1 = K^{11} \left( \frac{h_1 - h_2}{a} \right) + K^{12} \left( \frac{h_0 - h_1}{c} \right) + K^{13} \left( \frac{h_4 - h_1}{d} \right) \quad (127)$$

Similarly, the second component of the upscaled Darcy velocity (along the y direction) must be equal to the mean of the total flux passing through the composite block in that direction.

$$-V_2 = \frac{ac v_2^1 + bc v_2^2}{(a+b)e} \quad (128)$$

where $v_2^1$ represents the local scale Darcy velocity component in the second direction (the y direction) in the first block, and $v_2^2$ represents the local scale Darcy velocity component along the second direction for the second block. Replacing the local scale velocities by their expression, we obtain

$$-V_2 = \left[ \frac{b k_1^{11} k_2^{12} + a k_2^{11} k_1^{12}}{b k_1^{11} + a k_2^{11}} \right] \left( \frac{h_1 - h_2}{a} \right) + \frac{1}{a+b} \left[ ak_1^{22} + bk_2^{22} - \frac{ab(k_1^{12} - k_2^{12})^2}{b k_1^{11} + a k_2^{11}} \right] \left( \frac{h_0 - h_1}{c} \right) \\ + \frac{1}{a+b} \left[ ak_1^{23} + bk_2^{23} - \frac{ab(k_1^{12} - k_2^{12})(k_1^{13} - k_2^{13})}{b k_1^{11} + a k_2^{11}} \right] \left( \frac{h_4 - h_1}{d} \right) \quad (129)$$

A similar computation for the third component of the upscaled Darcy velocity gives:

$$-V_3 = \frac{ad v_3^1 + bd v_3^2}{(a+b)d} \quad (130)$$
to be compared with:

$$ -V_3 = \left[ \frac{b k_1^{11} k_2^{13} + a k_2^{11} k_1^{13}}{bk_1^{11} + ak_2^{11}} \right] \left( \frac{h_1 - h_2}{a} \right) \\ + \frac{1}{a+b} \left[ ak_1^{23} + bk_2^{23} - \frac{ab(k_1^{12} - k_2^{12})(k_1^{13} - k_2^{13})}{bk_1^{11} + ak_2^{11}} \right] \left( \frac{h_0 - h_1}{c} \right) \\ + \frac{1}{a+b} \left[ ak_1^{33} + bk_2^{33} - \frac{ab(k_1^{13} - k_2^{13})^2}{bk_1^{11} + ak_2^{11}} \right] \left( \frac{h_4 - h_1}{d} \right) \quad (131) $$

Finally, using these expressions, one can estimate all the components of the equivalent conductivity tensor as follows:

$$ K^{11} = \frac{(a + b) k_1^{11} k_2^{11}}{bk_1^{11} + ak_2^{11}} \quad (132) $$

$$ K^{22} = \frac{1}{a+b} \left[ ak_1^{22} + bk_2^{22} - ab \frac{(k_1^{12} - k_2^{12})^2}{bk_1^{11} + ak_2^{11}} \right] $$

$$ K^{33} = \frac{1}{a+b} \left[ ak_1^{33} + bk_2^{33} - ab \frac{(k_1^{13} - k_2^{13})^2}{bk_1^{11} + ak_2^{11}} \right] $$

$$ K^{12} = \frac{bk_1^{11}k_2^{12} + ak_2^{11}k_1^{12}}{bk_1^{11} + ak_2^{11}} $$

$$ K^{13} = \frac{bk_1^{11}k_2^{13} + ak_2^{11}k_1^{13}}{bk_1^{11} + ak_2^{11}} $$

$$ K^{23} = \frac{1}{a+b} \left[ ak_1^{23} + bk_2^{23} - ab \frac{(k_1^{12} - k_2^{12})(k_1^{13} - k_2^{13})}{bk_1^{11} + ak_2^{11}} \right] $$

These relations generalize the classical arithmetic and harmonic mean for a layered medium. One can check that Equation (132) reduces to a diagonal tensor with the arithmetic mean in directions y and z when the cross terms of the local conductivity tensors are zero. We also note that Equation (132) reduces to the equations that were presented by Kasap and Lake [105] for an analogous tensorial problem in 2D.

## 6. Numerical Experiments

Among the different techniques that were presented in the previous section on criteria for the equivalent conductivity, it is not obvious to know in advance which upscaling techniques (if any) might satisfy some other criteria in addition to those directly used for defining the equivalent conductivity.

For example, when one uses the DIAG method (criterion) with linearly varying head conditions (boundary conditions), a question that can be asked is the following: “is the TDP criterion also verified?”. Similarly, what about the Net Surface Flux when defining K using the Volume Average Flux? Other important questions are whether the resulting equivalent conductivity is a true second rank tensor, whether it is symmetric, and whether it is positive-definite? (Briefly, K is a “true” second rank tensor if its components $K_{ij}$ are transformed according to tensorial transformation rules under a rotation of the coordinate system; K is a non-symmetric tensor if $K^T \neq K$; and K is positive definite if its symmetric part $K^S = (K^T + K)/2$ is positive-definite. Finally, a symmetric matrix A is positive-definite if any quadratic form $x^T A x$ is strictly positive for all non-zero vectors x.)

In this section, we apply the different upscaling methods on numerical examples, and we explore the properties of the methods. The cases are chosen to illustrate not only the standard situations, but also the pathological situations that can occur, and to understand them. This series of numerical experiments provides counter examples allowing to know when some criteria or properties are not systematically satisfied. (In the next Section 7, we
will concentrate on analytical proofs to demonstrate the general validity of some numerical observations made in this section concerning the properties of the various versions of equivalent conductivity.)

## 6.1. Continuous and Composite Media in 2D

In the case of continuous and composite media, the detailed and global flow calculations, required to numerically compute the equivalent permeability tensor and to check the accuracy of the methods, were conducted either with the commercial finite element code FEFLOW (http://www.feflow.info, accessed in 18 June 2022) or with MATLAB routines developed by us for the purpose of this work. In both cases, we used standard bilinear rectangular finite elements.

### 6.1.1. Experimental Setup, Types of Media, and Methodology

Numerical tests have been conducted in two and three dimensions, but for illustration purposes we present only two-dimensional results. Here, we study 12 synthetic media with different spatial distributions of hydraulic conductivities. These media are represented in Figure 13.

<figure><img src="image_16.png" /><figcaption>**Figure 13.** The twelve two-dimensional synthetic samples of porous media proposed for analysis. The white and light grays represent high conductivity values, while the black and dark grays represent low conductivities.</figcaption></figure>

The first test sample, denoted K1, is a homogeneous anisotropic medium; the local conductivity $k_{ij}$ is a constant tensor represented by a matrix with non-zero of off-diagonal terms. All the other samples K2 to K12 are heterogeneous, and their local conductivity is isotropic, represented by a spatially varying scalar $k(x, y)$. However, in spite of local isotropy, some of these samples can exhibit anisotropic geometric structure due to the spatial organization of the local conductivities. Media K2 to K6 are stratified (or multilayered), with a contrast of conductivities of 100. Media K7 and K8 are multi-Gaussian random media, with a lognormal probability distribution of hydraulic conductivities, respectively,
with isotropic structure (K7) and anisotropic stratified structure (K8). The random medium K8 exhibits an inclined geometrical anisotropy (or statistical anisotropy). Media K9 and K10 are, respectively, isotropic and anisotropic random binary media, having the same proportions of the two phases. (The term “phases” is understood here, for binary media, in the sense of material phases (e.g., high permeability sand for the white phase and low permeability clay for the black phase). In all cases, the porous media are considered to be water-saturated (or saturated with a single fluid phase).) Medium K11 is a random Boolean medium made up of porous matrix (black) with “anisotropic” lenticular inclusions of elliptic shape. Medium K12 is similar to K6, but here, the unique conductive stratum connects two adjacent perpendicular faces (instead of two opposite parallel faces in K6).

For all these block samples of porous media, a reference hydraulic conductivity tensor can be estimated either with an exact or an approximate analytical formula (see Table 3). In all the cases, the reference hydraulic conductivity assumes that the medium is infinite and that the flow is uniform through the medium. For the stratified media, the reference equivalent conductivity is calculated using the arithmetic and harmonic means for the directions parallel and perpendicular to the strata.

<table><caption>**Table 3.** Reference conductivity tensors for the twelve two-dimensional synthetic media.</caption><thead><tr><th>Medium</th><th>Description</th><th>K<sub>11</sub> [m/s]</th><th>K<sub>22</sub> [m/s]</th><th>Angle [°]</th></tr></thead><tbody><tr><td>K<sub>1</sub></td><td>Homogeneous anisotropic</td><td>1.00 × 10<sup>-2</sup></td><td>1.00 × 10<sup>-4</sup></td><td>25</td></tr><tr><td>K<sub>2</sub></td><td>Stratified: 4 horizontal strata</td><td>5.05 × 10<sup>-3</sup></td><td>1.98 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>3</sub></td><td>Stratified: 8 horizontal strata</td><td>5.05 × 10<sup>-3</sup></td><td>1.98 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>4</sub></td><td>Stratified: 8 inclined strata</td><td>5.05 × 10<sup>-3</sup></td><td>1.98 × 10<sup>-4</sup></td><td>-45</td></tr><tr><td>K<sub>5</sub></td><td>Stratified: 16 inclined strata</td><td>5.05 × 10<sup>-3</sup></td><td>1.98 × 10<sup>-3</sup></td><td>-45</td></tr><tr><td>K<sub>6</sub></td><td>Single layer: 1 inclined stratum (layer)</td><td>9.12 × 10<sup>-3</sup></td><td>1.09 × 10<sup>-4</sup></td><td>-45</td></tr><tr><td>K<sub>7</sub></td><td>Random Multi-Gaussian isotropic</td><td>1.83 × 10<sup>-2</sup></td><td>1.83 × 10<sup>-2</sup></td><td>0</td></tr><tr><td>K<sub>8</sub></td><td>Random Multi-Gaussian anisotropic</td><td>2.48 × 10<sup>-2</sup></td><td>1.29 × 10<sup>-2</sup></td><td>25</td></tr><tr><td>K<sub>9</sub></td><td>Random Binary isotropic</td><td>9.33 × 10<sup>-4</sup></td><td>9.33 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>10</sub></td><td>Random Binary anisotropic</td><td>2.09 × 10<sup>-3</sup></td><td>4.49 × 10<sup>-4</sup></td><td>25</td></tr><tr><td>K<sub>11</sub></td><td>Random Binary Inclined ellipsoids</td><td>6.80 × 10<sup>-4</sup></td><td>1.70 × 10<sup>-4</sup></td><td>25</td></tr><tr><td>K<sub>12</sub></td><td>Single layer: 1 inclined stratum (layer)</td><td>8.30 × 10<sup>-3</sup></td><td>1.10 × 10<sup>-4</sup></td><td>35</td></tr></tbody></table>

For the random multi-Gaussian medium and for the binary *isotropic* medium with identical proportions of high and low permeability cells (K7 and K9, respectively), the exact effective conductivity is the geometric mean $K_G$ [30].

For the multi-Gaussian and binary *anisotropic* media (K8 and K10, resp.), we use the geometric mean corrected by the anisotropy factors proposed by Ababou [68] and Ababou ([106], Appendix B therein):

$$K_{ii} = K_G \exp \left[ \frac{\sigma_Y^2}{2} \left( 1 - \frac{2 l_h}{D l_i} \right) \right] \quad (133)$$

where $K_G$ represents the geometric mean, $\sigma_Y^2$ the variance of the neperian logarithm of conductivity, $l_i$ the correlation scale in the $i$ direction, and $l_h$ the arithmetic mean of the correlation scales. The same formula was also applied to the anisotropic Boolean medium K11 with elliptical inclusions. Finally, note that all media except K1 have *locally isotropic* conductivity $k(x, y)$; on the other hand, several of the proposed sample media such as K8, K10, and K11 are *structurally* or *statistically* anisotropic, although their local conductivity $k(x, y)$ is isotropic (scalar).

For each of these 12 samples of porous media, 8 equivalent conductivity tensors were computed numerically using the methods described in the previous section. These eight conductivity tensors correspond to four averaging methods—or interpretation methods
(DIAG, NSF, VSF, and VAF)—and to two types of boundary conditions (linearly varying head and permeameter type).

To analyze the results, we use several measures. First, we define a dimensionless error norm $\epsilon$ by computing the *Frobenius norm* of the difference between the estimated equivalent conductivity tensor and the reference tensor. The Frobenius norm is also called *Euclidian norm* [107]. The Frobenius or *Euclidian norm* was chosen because it can be used to compare two tensors globally. It is close to zero when the eigenvalues and eigenvectors of the two tensors are close. It is defined as follows:

$$ \|\mathbf{K}\|_{\text{FRO}} = \sqrt{\text{trace}(\mathbf{K}\mathbf{K}^{\text{T}})} \qquad (134) $$

where N = 2 in 2D, or N = 3 in 3D. The Frobenius norm can also be expressed as follows in terms of the eigenvalues of K, i.e., the principal conductivities $K_{11}^*, K_{22}^*, K_{33}^*$ (which are independent of the chosen coordinate system): $\|\mathbf{K}\|_{\text{FRO}} = \sqrt{K_{11}^{*2} + K_{22}^{*2} + K_{33}^{*2}}$. Interestingly, it can also be proved that the Frobenius norm of a symmetric definite positive tensor $K_{ij}$ is the trace of the “square-root matrix” B defined by $\mathbf{K} = \mathbf{B}\mathbf{B}^{\text{T}}$. In other words, we have the identity: $\|\mathbf{K}\|_{\text{FRO}} = \text{trace}(\mathbf{B})$. The proof is obtained by expressing B in the form $\mathbf{B} = \mathbf{L}\mathbf{D}^{1/2}$ where L is lower triangular and D is diagonal, and using the principal basis where $K_{ij}$ is represented by the diagonal matrix $K^*$.

Secondly, we define the error $\epsilon$ as the Frobenius norm of the difference between the “estimated” and the “reference” equivalent conductivity, and we normalize it by the norm of the reference conductivity, namely:

$$ \epsilon = \frac{\|\mathbf{K}_{\text{EST}} - \mathbf{K}_{\text{REF}}\|_{\text{FRO}}}{\|\mathbf{K}_{\text{REF}}\|_{\text{FRO}}} \qquad (135) $$

Thirdly, we estimate the relative magnitude $\beta$ of the anti-symmetry of the resulting equivalent tensor:

$$ \beta = \frac{\frac{1}{2} |K_{12} - K_{21}|}{\|K_{\text{FRO}}\|} \qquad (136) $$

Both criteria $\epsilon$ and $\beta$ are useful.

On the other hand, we use other comparison methods, based on the upscaled flow field (V); the estimated equivalent K tensor is used to compute again the flow field with the same boundary conditions that were used to compute the small-scale flow field (v). We also compare the total dissipated power and the fluxes obtained with the equivalent medium (equivalent conductivity K) and with the original heterogeneous medium (local conductivity k).

### 6.1.2. Comparisons with the Reference Permeability Tensors for 12 Porous Samples

In Figure 14, we first show for each 2D porous sample (K1, K2, ..., K12) the relative difference $\epsilon$ between the reference value and the computed tensor for all the techniques that we tested.

The first striking characteristic that can be seen in this graph is that there is no obvious and systematic classification of the methods. The volume surface flux with linearly varying head (VSF-GRD) seems to provide very often a tensor far from the reference values, but at the same time it performs better than many other techniques for the first case (medium K1). Some techniques seem to be more reliable than others, but the classification is often changing depending on the type of heterogeneous medium. None of the techniques systematically provides a value equal or very close to the reference.

It is therefore important to understand more precisely what are the physical mechanisms that control those differences. Let us start with the homogeneous anisotropic case (Test case 1 denoted “K1” in Table 3).
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Comparison between the twelve estimated ($K_{est}$) and reference ($K_{ref}$) hydraulic conductivity tensors for the 2D cases; the ordinate represents the normalized Frobenius norm of their difference, and the abscissa represents the different types of media (12 porous samples labeled by their identification number).</figcaption></figure>

* For the homogeneous anisotropic medium, when applying *linearly varying head boundary conditions*, all the methods correctly identify the permeability tensor.

* With permeameter (PRM) boundary conditions, only the VAF+VAG method provides the right answer (the reference flow field and the equivalent conductivity tensors are displayed in Figure 15a-c). We did not plot the DIAG method here since it is obvious that it will not properly capture the anisotropy. Since the homogeneous domain is anisotropic, the flow lines are not perpendicular to the head isovalues. The geometry of the flow is distorted in such a way that it minimizes the total dissipated power within the media accounting both for the boundary conditions and the anisotropy imposed by the conductivity tensor. As a consequence, there is a difference between the mean head computed on the two no-flow boundary conditions, and there is a deviated component of the mean Darcy velocities within the domain.

* When we compute the net surface flux (NSF), these inclined components cannot be identified since the two faces in the direction perpendicular to the vertical component are forced to be impervious. The head gradient in the domain is still correctly estimated, using the volume average gradient (VAG), and allows identifying an off-diagonal term in the equivalent conductivity tensor. As a result, $K_{NSF}$ correctly captures the direction of anisotropy (see Figure 15c), but the magnitude of the maximum principal value is underestimated and the perpendicular component is overestimated (the ellipsoid is closer to a circle than the reference).

* In contrast, the vectorial surface flux (VSF) averages directly the specific discharge vectors "v" on the faces of the sample and overestimates the vertical component because "v" tends to be very inclined along the boundary. Indeed, the "v" vector inside the domain (near no flow boundaries) is much less inclined but not sampled by the VSF. We observe that the VSF technique properly identifies the main directions of anisotropy, but in this case, it overestimates the maximum principal component and
underestimate the smallest component (the ellipsoid is flattened as compared to the reference one).

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Illustration of the head field (a) and velocity field (b) in a homogeneous anisotropic medium (sample K1) computed with permeameter type boundary conditions (PRM). The medium is homogeneous, as $K_1$, but it has been chosen here with a different orientation to better illustrate the distortion of the flow field. Its main axis of anisotropy in this figure is inclined at –45 degrees, $K_{max} = 10^{-4}$, $K_{min} = 10^{-5}$. Anisotropy ellipses are shown below in plot (c). These ellipses represent the square root of flux directional permeability $\sqrt{K_{flux}}$ [77] obtained here with the NSF, VSF, and VAF criteria. The ellipses axes are aligned with the principal directions of anisotropy, and their principal radii are the square roots of the principal values of conductivity.</figcaption></figure>

A brief conclusion from this first experiment (K1: homogeneous anisotropic medium) is the following: linearly varying head boundary conditions (GRD) and volume average flux (VAF) are best. However, this is not always true, depending on the type of medium, as will be seen in the next example.

We now examine the stratified or layered medium K2 displayed earlier in Figure 13. The results are shown in Figure 16. One can observe a behavior that is opposite to the one described in the previous paragraph. For test K2, all the techniques based on PRM boundary conditions are converging toward the reference value, whereas the techniques using GRD boundary conditions leads to a significant error. More precisely, for sample K2:

* all the techniques identify properly the principal directions of anisotropy;

* all the techniques also estimate correctly the largest principal component $K^{‖}$ parallel to strata (arithmetic mean of local permeabilities);

* in the perpendicular direction, we expect to obtain for $K^{⊥}$ the harmonic mean permeability; however, all the methods using GRD boundary conditions overestimate $K^{⊥}$ significantly (Figure 14); the overestimation is maximum for VSF, minimum for VAF, and NSF lies in between.
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Medium K2 with linearly varying head boundary conditions (GRD): (a) K field, the red layers are the high permeability layers. (b,c) Head distribution in the domain for the two flow experiments. (d) Ellipsoids of the square root of flux directional permeability obtained with the VAF, NSF, and VSF criteria. (e,f) Darcy velocity fields for the two flow experiments; the background color represents the norm of the Darcy velocity.</figcaption></figure>

The reason for this overestimation of $K^{\perp}$ is the presence of flow cells on the sides of the most permeable layers in the K2 sample (see Figure 16e,f). Water is entering the high permeability layers at their base and leaving them laterally because of the prescribed head gradient along the boundary. The same behavior would occur if the block would be immersed in a medium with intermediate conductivity: the flow would converge toward the high conductivity layer and diverge around the low conductivity zones. The consequence is that the total flux perpendicular to layers is much higher than it would be if the layers were infinite, and therefore we obtain a systematic overestimation of the equivalent conductivity $K^{\perp}$ perpendicular to layers.

Note: the opposite situation arises when forcing the flow through flux conditions instead of head conditions; this is discussed in Section 6.1.5 comparing flux immersion vs. head immersion conditions.

Is this apparent overestimation of the orthogonal $K^{\perp}$ really a problem? In fact, one could argue that it is not the block permeability $K^{\perp}$ that is overestimated, but rather that the reference used for comparison underestimates the true perpendicular block conductivity. Indeed, the reference used here assumes an infinite extension of layers, and/or a set of infinitely thin layers. In fact, the K2 test contains just a few thick layers of finite extension, and the local 2D configuration of the flow near the boundaries suggests that indeed, the equivalent perpendicular $K^{\perp}$ is truly larger than the harmonic mean in this particular case (Figure 16f). Overall, this discussion brings about the need for using skin immersion (or border regions) when computing the equivalent conductivity [55], e.g., in the case of thick layers with respect to block size. The magnitude of these differences will probably decrease when the number of layers increases (compare samples K2 and K3 displayed earlier in Figure 13); this will be due to a decrease in the size of the rotational flow cells, such as those observable at the intersection of domain boundary and layer interfaces in Figure 16f.

Continuing with the comparison of the computed permeability tensors versus reference conductivities, we now consider media K3 and K4 (which have inclined layers). We find as expected that the errors are maximal with the DIAG method compared with all other methods.
For the random multi-Gaussian fields (concerning the spatially correlated MultiGaussian random fields K7 and K8, the notions of *local* vs. *structural* or *statistical* anisotropy should be distinguished. Both K7 and K8 represent spatial distributions of *locally isotropic* conductivity *k*(*x*, *y*). However, the first one (K7) is also isotropic *structurally (statistically)*, while the second one (K8) is anisotropic *structurally (statistically)*.), the isotropic field (K7), and the anisotropic field (K8), which have fairly small spatial correlation lengths with respect to block size, all the techniques are reasonably good except the VSF+GRD method (i.e., Vectorial Surface Flux “VSF” combined with linearly varying head or fixed gradient condition “GRD”).

Finally, the highest differences between the computed equivalent conductivity and the reference values are observed for binary random media (K9, K10, K11) when using the VSF+GRD method (vectorial surface flux with linearly varying head conditions).

6.1.3. Analysis of the Degree of Anti-Symmetry of the Computed Tensors

Figure 17 shows for each numerical experiment (test media K₁, ..., K₁₂) the relative magnitude *β* of the anti-symmetric component of the resulting permeability tensors, normalized by the Frobenius norm of the full tensor (Equation (134)). It can be seen that anti-symmetry arises with some of the methods (not all), and only for some types of heterogeneity (not all).

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Magnitude *β* of the antisymmetric component of the equivalent permeability tensors, relative to the Frobenius norm of the full permeability tensor.</figcaption></figure>

We observe that anti-symmetry arises in the simplest test case (K1: homogeneous anisotropic medium) when applying the NSF method with permeameter (PRM) boundary conditions. This is not surprising, since the NSF method cannot deliver in a single experiment the two components of the average Darcy velocity when permeameter (PRM) boundary conditions are used.
For all the stratified media, horizontal or inclined (K2 to K5), all the methods gave symmetrical tensors. Those results are related to the geometric symmetry of the media themselves, with respect to heterogeneous structure as well as domain shape.

Indeed, non-symmetry arises for medium K6 containing a single inclined stratum that is not well centered into the block. In that situation, the only method that provides a symmetric tensor is the volume average flux VAF method with linearly varying head conditions, but all other methods yield a non-symmetric tensor.

Similarly, for the multi-Gaussian and binary media (K7, K8, K9, K10, K11), all the techniques except VAF produce non-symmetric permeability tensors.

To provide a synthesis of all the previous results obtained for various 2D media, we plot in Figure 18 the different upscaling techniques on a map in terms of the degree of non-symmetry. Each technique is represented by a point whose coordinates are as follows: the *x*-axis represents the average normalized error over all numerical experiments (twelve different media), and the *y*-axis represents the average of the normalized anti-symmetric component of the computed tensors. We see clearly in Figure 18 that the Volume Average Flux technique with the Linearly Varying head is the technique that always provides a symmetric tensor and that it also has the lowest error with respect to the reference value on average over all test cases (over the limited sample of 12 blocks of porous media investigated here).

<figure><img src="image_21.png" /><figcaption>**Figure 18.** Synthetic representation of the comparison between the reference and the numerically estimated equivalent conductivity tensors, in terms of the degree of anti-symmetry. The abscissa of a point in this figure represents the average of the norm over the twelve 2D blocks of porous media defined in the text. The ordinate represents the mean degree of anti-symmetry, averaged over the twelve 2D media.</figcaption></figure>

### 6.1.4. Net Surface Flux and Total Dissipated Power Criteria (Verification)

In this section, we check numerically whether the Net Surface Flux and the Total Dissipated Power criteria are verified a *posteriori*. For that purpose, we solve the flow problem numerically on the upscaled media. Since some approximations had to be made in order to define the equivalent *K* with some techniques (NSF for example), it is not obvious that the criteria are verified a *posteriori*.
In Figure 19, we compare the net surface flux in the $x_1$ direction computed on the detailed heterogeneous media versus the corresponding equivalent homogeneous media with the PRM and GRD boundary conditions, and this was performed for all the media (the 12 porous block samples). The following results emerge from these graphical comparisons:

* The DIAG technique reproduces very well the Net Surface Flux in the direction of flow imposed by the boundary condition. Note, in the direction perpendicular to the main gradient (GRD case), the fluxes are necessarily zero because the DIAG method neglects the off-diagonal term of the conductivity tensor.

* The NSF method is almost perfect, but still one can see that with PRM boundary conditions there is a slight discrepancy due to the approximation discussed earlier.

* The VSF method shows the highest dispersion and seems to be the least accurate.

* The VAF method seems good overall but shows some slight difference between the fluxes in the case of the GRD boundary conditions. Surprisingly, it seems that the VAF perfectly reproduces the fluxes in the case of PRM conditions.

<figure><img src="image_22.png" /><figcaption>**Figure 19.** Comparison between the values of the microscopic total flux $Q_1$ calculated on the heterogeneous media (abscissa) with the total macroscopic flux computed with the same boundary conditions on the upscaled homogeneous medium (ordinate). Each graph represents all the numerical results for a given upscaling method (DIAG, NSF, VSF, VAF).</figcaption></figure>

Similar comparisons are also shown in terms of Total Dissipated Power (TDP) in Figure 20. It can be seen that, among all the techniques that we tested, only the VAF method accurately reproduces the TDP criterion. Among the other techniques, the DIAG and NSF methods perform rather well overall, while the VSF method is less accurate.
<figure><img src="image_23.png" /><figcaption>**Figure 20.** Comparison between the total microscopic dissipated power (P) computed on the heterogeneous media (abscissa) with the total macroscopic power calculated with same boundary conditions on the homogeneous media (ordinate). Only the VAF method with PRM and GRD boundary conditions ensures the equality of the dissipated powers.</figcaption></figure>

All these numerical/graphical results are also summarized in Table 4.

### 6.1.5. Flux Immersion Conditions: Comparisons with Linear Head Immersion

This section discusses prescribed “flux immersion” boundary conditions versus the prescribed head conditions. More precisely, “flux immersion” conditions correspond to the case where the flow is forced by immersion of the porous block in a flux vector field (using the normal fluxes as BC’s)—to be compared to hydraulic gradient immersion (where linearly varying heads are imposed on the boundary).

We use a horizontally layered block K2bis as a prototype case for this discussion. This block is similar to a thickly fractured block comprising three layers (M; F; M), with the “Fracture” layer sandwiched between the “Matrix” layers as illustrated earlier in Figure 11. The relative permeabilities and thicknesses of the three layers are ($k_M$; $k_F$; $k_M$) = (1; 1000; 1) and ($e_M$; $e_F$; $e_M$) = (9; 2; 9). The averaging method used to interpret the flux immersion experiment is VAF for the flux and VAG for the gradient. The results obtained with flux immersion $K_{FLX}$, compared to linear head or “gradient immersion” $K_{GRD}$, are the following:

$$ K_{FLX} \approx \begin{bmatrix} 2.3 & 0 \\ 0 & 1.11 \end{bmatrix} \qquad (137) $$

$$ K_{GRD} \approx \begin{bmatrix} 100.9 & 0 \\ 0 & 10.8 \end{bmatrix} \qquad (138) $$
<table><caption>**Table 4.** Main properties of the upscaling techniques discussed in this paper.</caption><thead><tr><th rowspan="2">Upscaling Method</th><th rowspan="2">Boundary conditions</th><th colspan="6">Properties of the Method</th></tr><tr><th>NSF</th><th>TDP</th><th>Symmetry</th><th>Reproduction of Reference</th></tr></thead><tbody><tr><td rowspan="2">Averaging techniques</td><td rowspan="2">PRM</td><td>Analytical proof</td><td>Numerical test</td><td>Analytical proof</td><td>Numerical test</td><td>Analytical proof</td><td>Numerical test</td></tr><tr><td>BD</td><td>Smiley</td><td>Smiley</td><td>Smiley</td><td>Diagonal tensor</td><td>Smiley</td></tr><tr><td colspan="2"></td><td colspan="6">§7.1</td></tr><tr><td>DIAG</td><td>-</td><td>GRD</td><td>NN</td><td>Smiley</td><td>Smiley</td><td>Diagonal tensor</td><td>Smiley</td></tr><tr><td>NSF</td><td>VAG</td><td>PRM</td><td>NN</td><td>Smiley</td><td>Smiley</td><td>NN</td><td>Smiley</td></tr><tr><td>NSF</td><td>-</td><td>GRD</td><td>BD</td><td>Smiley</td><td>Smiley</td><td>NN</td><td>Smiley</td></tr><tr><td>VSF</td><td>VAG</td><td>PRM</td><td>NN</td><td>Smiley</td><td>Smiley</td><td>NN</td><td>Smiley</td></tr><tr><td>VSF</td><td>-</td><td>GRD</td><td>NN</td><td>Smiley</td><td>Smiley</td><td>NN</td><td>Smiley</td></tr><tr><td>VAF</td><td>VAG</td><td>PRM</td><td>Smiley</td><td>Smiley</td><td>Smiley</td><td>NN</td><td>Smiley</td></tr><tr><td>VAF</td><td>-</td><td>GRD</td><td>NN</td><td>Smiley</td><td>Smiley</td><td>7.2</td><td>Smiley</td></tr><tr><td>-</td><td>VAG</td><td>FLX</td><td>BD</td><td>Smiley</td><td>Smiley</td><td>Smiley</td><td>Smiley</td></tr><tr><td>VAF</td><td>VAG</td><td>SKN</td><td colspan="5">NN</td></tr><tr><td colspan="2"></td><td colspan="5">[55]</td></tr><tr><td>VAF</td><td>-</td><td>PRD</td><td>BD</td><td>Smiley</td><td>Smiley</td><td>Smiley</td><td>Smiley</td></tr><tr><td colspan="2"></td><td colspan="5">[96]</td></tr><tr><td colspan="2"></td><td colspan="5">[96]</td></tr><tr><td colspan="2"></td><td colspan="5">[96]</td></tr><tr><td colspan="2"></td><td colspan="5">[37]</td></tr></tbody></table>

**Legend:** ☺ = The criterion is always verified. ☐ = The criterion is not always verified. ⊗ = The criterion is strongly not respected. NR = Not Relevant (or criterion cannot be computed). NN = Not Necessary (a sufficient counter example is provided). BD = The criterion is imposed By Definition.

These results can be interpreted by comparing them with the reference tensor:

$$ \mathbf{K}_{REF} = \begin{bmatrix} \bar{K}_A & 0 \\ 0 & \bar{K}_H \end{bmatrix} \approx \begin{bmatrix} 100.9 & 0 \\ 0 & 1.11 \end{bmatrix} \quad (139) $$

Note that the equivalent block permeability is expressed in the reference frame aligned with the layers, with $K_A$ and $K_H$ denoting the weighted arithmetic and harmonic means of the permeabilities of the three layers (The arithmetic mean “$K_A$” should not be confused with the anti-symmetric part also denoted “$K_A$” in other sections.). This reference tensor was shown in Section 5.1 to be the exact equivalent permeability of the block under piece-wise linear head conditions, or equivalently, under piecewise constant gradient conditions (and therefore also under piecewise constant flux conditions as well). Comparisons with the GRD and FLX numerical experiments show that:

* gradient immersion (GRD) significantly overestimates the permeability orthogonal to layers ($K^{\perp} \approx 10.8 \gg 1.11$);
■ flux immersion (FLX) drastically underestimates the permeability parallel to layers ($K'' \approx 2.3 \ll 100.9$).

These results are illustrated in Figure 21 and they confirm experimentally that the two tensors are ordered by the relation $K_{FLX} \le K_{GRD}$ in the following sense: matrix $A \le B$ if matrix $B - A \ge 0$ (B-A is semi-definite positive). Pouya et al. [43,51,85] proved analytically in the continuous case the same inequality $K_{FLX} \le K_{GRD}$.

<figure><img src="image_24.png" /><figcaption>**Figure 21.** Schematic representation of the 3 permeability tensors ellipsoids for a finite 3 horizontal layers block (matrix/fracture block).</figcaption></figure>

Furthermore, the tensors $K_{GRD}$ (obtained under gradient immersion conditions, “conductive upscaling”) and $K_{FLX}$ (obtained under flux immersion conditions, “resistive” upscaling) are bounding the reference tensor $K_{REF}$.

From a theoretical point of view, the flux immersion and hydraulic gradient immersion experiments can be analyzed in the light of other theoretical approaches, the first one (A) dealing directly with a flux governing equation to analyze local flow, and the second one (B) considering together the conductivity upscaling problem (under gradient conditions) and the “dual” resistivity upscaling problem under flux conditions [75]. These approaches were outlined earlier in Sections 2.2.2 and 4.3 and are also described in detail in Appendix A.2.

However, in the light of the resistivity approach from Fadili and Ababou [75], the result $K_{FLX} \le K_{GRD}$ can be re-interpreted as $[Rij]^{-1} \le [Kij]$ if we choose to denote $[Rij]^{-1}$ the upscaled resistivity (under prescribed velocity) and Kij the upscaled conductivity (under prescribed gradient, or linearly varying head). Thus, if we admit that the equivalent [Rij] 1 and Kij tensors correspond to two distinct flow experiments with boundary conditions FLX and GRD (respectively), then, the two tensors do not necessarily match, and they are always ordered as indicated above.

## 6.2. Simple 2D Fracture Networks

In this section, we use simple 2D fracture networks (Figures 22 and 23) to illustrate further the properties of the different techniques. Only linearly varying head boundary conditions are used in this section.

<table><caption>**Table 5.** Hydraulic conductivity tensors computed on the 2D fracture networks displayed in Figures 22 and 23.</caption><thead><tr><th></th><th colspan="2">Case (a)</th><th colspan="2">Case (b)</th><th colspan="2">Case (c)</th><th colspan="2">Case (d)</th></tr></thead><tbody><tr><td>K<sub>NSF</sub></td><td>(1.42 0.84)<br>(0.84 1.42)</td><td>(1.42 0.84)<br>(0.0 1.42)</td><td>(1.42 0.84)<br>(0.0 1.42)</td><td>(1.42 0.84)<br>(-0.84 1.28)</td><td>(1.42 0.84)<br>(0.78 1.28)</td><td>(1.42 0.84)<br>(1.53 1.28)</td><td>(1.42 0.84)<br>(1.53 1.28)</td></tr><tr><td>K<sub>VSF</sub></td><td>(1.63 1.44)<br>(1.44 1.63)</td><td>(1.42 0.72)<br>(0.72 1.63)</td><td>(1.63 0.0)<br>(0.0 1.63)</td><td>(1.63 0.0)<br>(0.0 1.63)</td><td>(1.53 1.28)<br>(1.44 1.59)</td><td>(1.53 1.28)<br>(1.44 1.59)</td></tr><tr><td>K<sub>VAF</sub></td><td>(1.91 1.68)<br>(1.68 1.91)</td><td>(1.65 0.84)<br>(0.84 1.91)</td><td>(1.91 0.0)<br>(0.0 1.91)</td><td>(1.91 0.0)<br>(0.0 1.91)</td><td>(1.89 1.63)<br>(1.63 1.89)</td><td>(1.89 1.63)<br>(1.63 1.89)</td></tr></tbody></table>
<figure><img src="image_25.png" /><figcaption>**Figure 22.** Schematic of a fracture network and the Darcy velocities corresponding to a flow experiment with a prescribed horizontal gradient (**a**–**c**) or a prescribed vertical gradient (**d**–**f**). Figures (**a**, **d**) represent the network and the two sets of boundary conditions. Figures (**b**, **e**) represent the velocity vectors in the individual fractures. Figures (**e**, **f**) represent the prescribed mean gradient and the resulting mean velocity. In both cases, the Darcy velocity is deviated toward the principal direction of anisotropy, which is oriented at 45 degrees.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 23.** Schematic of three simple fracture networks (**b**–**d**), with the Darcy velocities corresponding to a flow experiment with a constant horizontal gradient. All the fractures have the same aperture. The lengths of the 4 fractures are identical in cases **b** and **c**, and the local conductivities are also identical. The letters (**b**–**d**) above the graphs correspond to the cases discussed in the text and to the results reported in Table 5.</figcaption></figure>

First, we want to illustrate when non symmetry occurs. For that purpose, we start
with a simple medium (Figure 22) such that all the techniques deliver a symmetric tensor
with a strong off-diagonal component (Table 5, case a). After having described this medium
and the corresponding results, we will modify it to show when the non-symmetry occurs.

The initial medium contains two fractures (Figure 22) intersecting exactly in the center of the domain. The medium is symmetric around a line going through the center of the domain and inclined at 45 degrees. On Figure 22a, we see the Darcy velocities computed for a linear gradient oriented along the x₁ direction. J is horizontal; the average velocity cannot be in the same direction due to the anisotropy and is deviated toward the principal direction of anisotropy at 45 degrees (Figure 22b). When the prescribed gradient is rotated by 90 degrees (Figure 22c), the velocity is again deviated toward the principal direction
of anisotropy at 45 degrees (Figure 22d). Due to the geometrical symmetry of the fracture system, the flux $\langle v_2 \rangle_{NSF}^1$ is identical to $\langle v_1 \rangle_{NSF}^2$, and this directly ensures the symmetry of $K_{NSF}$ (see Table 5, case a). Still because of the geometrical symmetry, the vectorial surface fluxes $\langle v_2 \rangle_{VSF}^1$ and $\langle v_1 \rangle_{VSF}^2$ are identical; $K_{VSF}$ is symmetric too as well as $K_{VAF}$ (Table 5, case a).

We now consider the fracture network (b) on Figure 23. Compared to the previous figure (Figure 22), here one boundary point of the upper left fracture has been moved toward the left. This has been done to ensure that the total flux in the vertical direction $\langle v_2 \rangle_{NSF}^1$ becomes zero during the experiment with a horizontal gradient, while keeping a value $\langle v_1 \rangle_{NSF}^2$ different from zero in the perpendicular experiment. As expected, the result is that $K_{NSF}$ is asymmetric (Table 5, case b). The other averaging techniques still provide a symmetric tensor as shown in Table 5.

If we continue to perturb the initial setup and move now the lower right fracture in order to obtain a perfect cross with perpendicular fractures but not aligned with the coordinates' axis (Figure 23c), the asymmetry becomes even stronger for $K_{NSF}$ with opposite off diagonal terms 0.84 and -0.84 (Table 5, case c). These results are obtained because when the gradient is prescribed in the horizontal direction, the fractures are oriented at an angle close to 40 degrees upward but the net fluxes on the faces are such that the estimated mean velocity with the NSF method is oriented downward. When the gradient is rotated by 90 degrees, again the mean velocity estimated with the NSF criteria is inclined toward the right. In both conditions, the velocity field always rotates with respect to the prescribed gradient, and the only way to describe that behavior is to use an anti-symmetric tensor. We also note that in this example, the symmetric part of $K_{NSF}$ is perfectly diagonal. This is in close agreement with the other averaging techniques which deliver a diagonal tensor (Table 5) and therefore consider the medium as isotropic.

With the previous examples, $K_{VSF}$ and $K_{VAF}$ was always symmetric. However, in some situations, it happens that the tensor becomes asymmetric but with a moderate anti-symmetric component. To illustrate that behavior, we show an example in Figure 23 and Table 5 case (d); it is based on the initial case (a) displayed in Figure 22, except that the central node was moved in order to break the symmetry of the fracture network. In this new situation (Figure 23 and Table 5 case (d)), both the NSF and VSF criteria give a non-symmetric K tensor. Only $K_{VAF}$ is symmetric.

The next example (Figure 24a) illustrates when a negative principal value can occur with the VSF method. The fractured sample is symmetric around a horizontal line splitting the domain in two equal halves. In each half domain, a two-part fracture enters the domain from the left and leaves it from the side. When prescribing GRD boundary conditions, we impose a constant head difference between the entry and outlet of the fracture. The flux through the fracture is constant, but the orientations of the inflowing and outflowing vectors are different and even almost opposite. More precisely, the orientation of the flux entering the block from the left is more inclined than the orientation of the outflowing edge (Figure 24a). As a consequence, the sum of the two vectors gives a backward horizontal component (Figure 24b). Due to the symmetry, the sum of those vectors for the two fractures gives a zero vertical component. Overall, while the prescribed gradient should force a forward flux, the VSF method computes a mean flux that is oriented backward and therefore leads to a negative value of $K_{11}$ (Table 6).

<table><caption>**Table 6.** Hydraulic conductivity tensors computed on the networks displayed in Figure 24.</caption><thead><tr><th>K<sub>NSF</sub></th><th>K<sub>VSF</sub></th><th>K<sub>VAF</sub></th></tr></thead><tbody><tr><td>(9.7 0.0)<br>(0.0 7.3)</td><td>(-0.4 0.0)<br>(0.0 6.1)</td><td>(7.8 0.0)<br>(0.0 4.4)</td></tr></tbody></table>
<figure><img src="image_27.png" /><figcaption>**Figure 24.** Schematic of a fracture network leading to a negative diagonal component for VSF. (a) Geometry of the network and Darcy velocities for the two upper edges of the network when a horizontal gradient is prescribed. (b) Geometrical representation of the sum of the two vectors represented in subplot (a).</figcaption></figure>

On the other hand, in the same situation the forward flux is estimated properly by the Net Surface Flux and the Volume Average Flux methods, and both methods yield a positive definite K, with more or less similar positive values for $K_{11}$ (Table 6).

To conclude, this example shows that the VSF method does not necessarily provide a positive-definite tensor.

In the last example (Figure 25), we consider two set of fractures having a fixed spacing and fixed orientations.

<figure><img src="image_28.png" /><figcaption>**Figure 25.** Schematic of a fracture network made of two sets of equally spaced and parallel fractures.</figcaption></figure>

Based on this geometry, we analytically compute the reference equivalent permeability tensor for an infinite medium [59,60]. Briefly, this reference permeability is obtained by a superposition of fluxes (Q in m²/s). In 2D, the direct superposition of fluxes projected on the fracture segments yields, for each fracture f:

$$Q_f = \frac{g}{12\nu} b_f^3 \{ \mathbf{J} - (\mathbf{J} \cdot \mathbf{n}_f) \mathbf{n}_f \} \quad (140)$$

Summing over all fractures ($f = 1, \dots, N$) and using indicial notation, this yields:

$$Q_i = T_{ij} J_j = \frac{g}{12\nu} \sum_f b_f^3 \{ \delta_{ij} - (n_f)_i (n_f)_j \} J_j \quad (141)$$

The equivalent tensorial transmissivity $T_{ij}$ (m²/s) emerges). Finally, dividing by domain cross-sectional area (or domain length in the 2D case), e.g., for the case of a single
set of parallel fractures, there emerges a mean fracture spacing L, and one finally obtains (for the single set of parallel fractures):

$$q_i = K_{ij}I_j = \frac{g}{12\nu} \sum_f \frac{b_f^3}{L} \left\{ \delta_{ij} - (n_f)_i (n_f)_j \right\} J_j \quad (142)$$

where $K_{ij}$ is the equivalent hydraulic conductivity, obtained here explicitly. For a set of two parallel families of fractures, these formulas specialize as follows:

$$K = \frac{\rho g}{12\mu} \sum_{m=1}^{N} \frac{e_m^3}{L_m} (I - n_m n_m^T) \quad (143)$$

where N is the number of families of fractures, I is the identity matrix, for each family of fracture, $n_m$ is the unit vector normal to the fractures (they are all parallel), $n_m^T$ its transpose, $L_m$ is the spacing between two fractures, and $e_m$ the aperture.

By comparing the analytical results with the numerical ones (Table 7), we find that all the components of the tensor are correctly estimated only if we use the VAF method. The NSF method tends to underestimate the off-diagonal terms and the smallest principal value. The VSF method overestimates all the components of the tensor.

<table><caption>**Table 7.** Hydraulic conductivity tensors computed on the networks displayed in Figure 25.</caption><thead><tr><th>K<sub>NSF</sub></th><th>K<sub>VSF</sub></th><th>K<sub>VAF</sub></th><th>K<sub>ref</sub></th></tr></thead><tbody><tr><td>(1.4 0.5)</td><td>(1.5 0.7)</td><td>(1.4 0.6)</td><td>(1.4 0.6)</td></tr><tr><td>(0.5 0.5)</td><td>(0.7 0.7)</td><td>(0.6 0.6)</td><td>(0.6 0.6)</td></tr></tbody></table>

## 7. Analytical and Algebraic Proofs of the Properties of Equivalent K

In the previous section, we have seen some systematic patterns in our numerical results. While some techniques were shown that do not provide symmetrical tensors or do not satisfy the criterion of total dissipation of energy, others were always satisfying those criteria. These results are summarized in Table 5. In this section, we provide some analytical studies and proofs of the properties of the equivalent permeability tensors defined previously. The results of this work are also listed in Table 4 to confirm, when possible, the numerical observations.

### 7.1. DIAG + PRM ensures NSF and TDP

In the case of permeameter type boundary conditions, the Total Dissipated Power (TDP) criterion can be expressed very simply as follows (P is the dissipated power):

$$P = \rho g \Delta h_1 q_{B1} \quad (144)$$

This relation is derived in detail in the Appendix A.7.4 by replacing the volume integral in the definition of the TDP by a surface integral. All the terms on the no-flow boundaries vanish and only the terms on the boundaries with a prescribed head remain. Since the head is constant on those boundaries, it can be taken out of the integral and the whole equation simplifies to Equation (144) both for the heterogeneous and homogeneous media. The DIAG method ensures that the total flux in the direction of the imposed gradient is identical in the homogeneous and heterogeneous media $q_{B1} = Q_{B1}$, and the prescribed head difference $\Delta h_1$ between inlet and outlet is also identical on both media $\Delta h_1 = \Delta H_1$. Therefore, the total dissipated power criterion is always verified when we compute the equivalent conductivity using the DIAG + PRM method.

Of course, this technique does not allow computing the off-diagonal components of the tensor and therefore the result may be highly inaccurate if the principal direction of anisotropy of the heterogeneous medium are not aligned with the faces of the block. In the laboratory, the presence of anisotropy can be tested by measuring the heads on the sides of
the block and computing the head gradient in the direction perpendicular to the main flow. If this gradient is not equal to zero, then it is a clear indication that the sample has not been well oriented and that the resulting diagonal tensor will be inaccurate (Renard et al., 2001).

## 7.2. VAF + GRD Implies TDP

In the case of Linearly Varying Head boundary conditions, the total dissipated power can be expressed for any geometry as follows (see Table 2 and Appendix A.7.2):

$$P = -\rho g \mathbf{J} \cdot \langle \mathbf{v} \rangle_{\text{VAF}} \qquad (145)$$

As described in Renard [80] or Pouya and Fouche [43], the same formula can be applied both for the microscale and the homogenized macroscale problems; also, because the volume average flux criterion with linearly varying head enforces $\langle \mathbf{v} \rangle_{\text{VAF}} = \langle \mathbf{v} \rangle_{\text{VAF}}$, a similar equality holds for the dissipated power as a direct consequence of the VAF criterion:

$$P = -\rho g (\mathbf{J} \cdot \langle \mathbf{v} \rangle_{\text{VAF}}) = -\rho g (\mathbf{J} \cdot \langle \mathbf{V} \rangle_{\text{VAF}}) \qquad (146)$$

This last relation states the equality of the volume averaged dissipated power P of (i) the microscale medium and (ii) the homogenized macroscale medium.

## 7.3. Symmetry of K for VAF + GRD

Using VAF for averages (Volume Averaged Flux) and GRD for boundary conditions (gradient, i.e., linearly varying head condition), together, ensures the symmetry of the equivalent conductivity tensor K.

This was shown by several authors including Pouya and Fouche [43] for fractured media. We briefly present the proof in a slightly different form, with a tensorial local conductivity $k(\mathbf{x})$. The proof is based on the principle that a tensor K is symmetric if for any pair of vectors J and J', the scalar products $\mathbf{J} \cdot \mathbf{K} \mathbf{J}'$ and $\mathbf{J}' \cdot \mathbf{K} \mathbf{J}$ are equal. (Note that J' is just another vector, it is not the transpose of J. In this text, the notation for the transpose a vector J is $\mathbf{J}^T$ not $\mathbf{J}'$.)

To obtain the result, we consider two flow experiments in which J and J' are the prescribed boundary conditions, h and h' are the head solutions of the flow problems for those two boundary conditions, and v and v' are the corresponding Darcy velocities. The local hydraulic conductivity k is assumed to be a symmetric tensor, and therefore we can write the following equality:

$$\nabla h \cdot (\mathbf{k} \nabla h') = \nabla h' \cdot (\mathbf{k} \nabla h) \qquad (147)$$

By Darcy's law, the left term of the previous equation can be written as the divergence of the product $h\mathbf{v}$:

$$\nabla h \cdot \mathbf{k} \nabla h' == -\nabla h \cdot \mathbf{v}' = -\nabla (h \mathbf{v}') \qquad (148)$$

Now, since we are looking at the VAF method, we integrate the above equation over the domain of interest:

$$-\frac{1}{\Omega} \int_{\Omega} \nabla h \cdot (\mathbf{k} \nabla h') d\omega = \frac{1}{\Omega} \int_{\Omega} \nabla (h \mathbf{v}') d\omega = \frac{1}{\Omega} \int_{\Gamma} h (\mathbf{v}' \cdot \mathbf{n}) d\gamma \qquad (149)$$

Additionally, since we are using the GRD type of boundary condition $h(\mathbf{x}) = h_0 - \mathbf{J} \cdot \mathbf{x}$, we can insert it as follows (recall that the global hydraulic gradient J is defined as $\mathbf{J} = -\nabla h = -\nabla H$):

$$-\frac{1}{\Omega} \int_{\Omega} \nabla h \cdot (\mathbf{k} \nabla h') d\omega = \frac{1}{\Omega} \int_{\Gamma} (h_0 - \mathbf{J} \cdot \mathbf{x})(\mathbf{v}' \cdot \mathbf{n}) d\gamma \\ = h_0 \left( \frac{1}{\Omega} \int_{\Gamma} \mathbf{v}' \cdot \mathbf{n} d\gamma \right) - \mathbf{J} \cdot \left( \frac{1}{\Omega} \int_{\Gamma} \mathbf{x} (\mathbf{v}' \cdot \mathbf{n}) d\gamma \right) \qquad (150)$$