# Direct simulation of non-additive properties on unstructured grids

Pauline Mourlanette, Pierre Biver, Philippe Renard, Benoit Noetinger,
Guillaume Caumon, Yassine Alexandre Perrier

► To cite this version:

Pauline Mourlanette, Pierre Biver, Philippe Renard, Benoit Noetinger, Guillaume Caumon, et al. Direct simulation of non-additive properties on unstructured grids. Advances in Water Resources, 2020, 143, pp.103665. 10.1016/j.advwatres.2020.103665. hal-02891780

HAL Id: hal-02891780

https://hal.univ-lorraine.fr/hal-02891780

Submitted on 7 Jul 2020

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 Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Direct simulation of non-additive properties on unstructured grids

Pauline Mourlanette<sup>a,b,c,*</sup>, Pierre Biver<sup>a</sup>, Philippe Renard<sup>b</sup>, Benoît Nøtinger<sup>d</sup>, Guillaume Caumon<sup>e</sup>, Yassine Alexandre Perrier<sup>f</sup>

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

## ARTICLE INFO

**Keywords:**
* Upscaling
* Permeability
* Geostatistics
* Power-averaging

## ABSTRACT

Uncertainties related to permeability heterogeneity can be estimated using geostatistical simulation methods. Usually, these methods are applied on regular grids with cells of constant size, whereas unstructured grids are more flexible to honor geological structures and offer local refinements for fluid-flow simulations. However, cells of different sizes require to account for the support dependency of permeability statistics (support effect).

This paper presents a novel workflow based on the power averaging technique. The averaging exponent $\omega$ is estimated using a response surface calibrated from numerical upscaling experiments. Using spectral turning bands, permeability is simulated on points in each unstructured cell, and later averaged with a local value of $\omega$ that depends on the cell size and shape.

The method is illustrated on a synthetic case. The simulation of a tracer experiment is used to compare this novel geostatistical simulation method with a conventional approach based on a fine scale Cartesian grid. The results show the consistency of both the simulated permeability fields and the tracer breakthrough curves. The computational cost is much lower than the conventional approach based on a pressure-solver upscaling.

## 1. Introduction

Subsurface phenomena cannot be observed directly due to scale issues and inaccessibility. In hydrogeology, the high spatial variability of rock types and the associated permeability field as well as the high spatial and temporal variability of fluid types and displacements is a main source of uncertainties. Estimating these uncertainties is particularly important in a context of engineering design and decision making. Examples of applications include the management of over-exploited aquifers or the propagation of dangerous contaminants (De Marsily et al., 1998). Similar uncertainty issues coupled with decision-making are encountered in other applied geoscience engineering, such as oil and gas industry (Preux, 2016), CO₂ storage in aquifers (Michael et al., 2010; Akhurst et al., 2015), or geothermal energy production (Vogt et al., 2010; Quinlivan et al., 2015; Witter et al., 2019).

This paper focuses on the spatial modeling of permeability fields. Due to rock heterogeneity, a measurement at point-support is not sufficient to determine accurately the value at another point. To assess the

uncertainty related to this problem, one can employ geostatistical simulation methods (Deutsch and Journel, 1992; Goovaerts, 1997; Chilès and Delfiner, 2012).

To represent the spatial variability and to solve numerically the flow and transport equations, a mesh is used and its cells are populated with petrophysical properties. Classical meshes, widely known as regular structured (or Cartesian) grids, are composed of hexahedral cells in three dimensions. The interesting aspect of regular structured grids is that their cells all share the same size. Depending on the type of numerical methods employed to solve the flow problem, other types of grids exist. Some tools use irregular structured grids, also called corner-point (or stratigraphic) grids, obtained by distorting hexahedra to follow more faithfully geological layers. Regular and irregular structured grids are practical for the rapidity of calculation but they present some drawbacks: some cells can be degenerated for high thickness variations, the geometry of the intersections between faults or wells and the geological structure are approximated, and corner-point regular grids lack spatial adaptivity to honor local complex features of the parameter field. To address these issues, unstructured grids can be used: each vertex can

<footnote><sup>a</sup> Total, Avenue Larribau, 64000 Pau, France</footnote>

<footnote><sup>b</sup> University of Neuchâtel, Center for Hydrogeology and Geothermics, 11 Rue Emile Argand, CP-2000 Neuchâtel, Switzerland</footnote>

<footnote><sup>c</sup> Sorbonne University, 21 Rue de l'Ecole de médecine, 75006 Paris, France</footnote>

<footnote><sup>d</sup> IFPEN, 1-4 Avenue du Bois Préau, 92852 Rueil-Malmaison, France</footnote>

<footnote><sup>e</sup> Université de Lorraine, RING GeoRessources Laboratory, 2 rue du Doyen Marcel Roubault, Vandoeuvre-lès-Nancy, 54500, France</footnote>

<footnote><sup>f</sup> Université de Lorraine, ENSG, 2 rue du Doyen Marcel Roubault, Vandoeuvre-lès-Nancy, 54500, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: pauline.mourlanette@orange.fr (P. Mourlanette), pierre.biver@total.com (P. Biver), philippe.renard@unine.ch (P. Renard), benoit.noetinger@ifpen.fr (B. Nøtinger), guillaume.caumon@ensg.fr (G. Caumon).</footnote>
have a different number of neighbors, and the cell's sizes can vary. Indeed, the areas of the grid that present more interest or stronger variations for flow estimations can be refined and other parts coarsened (Prévost et al., 1996). Usually, to minimize numerical errors, the grid is refined around wells or faults but it can also be optimized based on flow patterns (Mlacnik et al., 2003). In addition, the cells can have various shapes such as tetrahedron or Voronoï polyhedron, offering more flexibility to adapt to geological heterogeneity (Blessent et al., 2011; Merland et al., 2014).

However, unstructured grids brings an additional difficulty for geostatistical simulations. Indeed, the statistics of the permeability (and of the other variables) are support dependent (i.e., in stationary settings, smaller cells tend to have a smaller internal variability and larger intercell variability). This support effect is well documented both theoretically and experimentally (Matheron, 1967; Dagan, 1993; Tidwell and Wilson, 1997): the probability distribution function of the permeability is different for cells of different volumes, as well as the spatial correlation. Intuitively, small cells will follow the geostatistical properties defined at the point-support scale, keeping their local average and covariance. In opposition, for very large blocks having a size greater than the underlying integral scale, the associated permeability value may be shown to stabilize to the so-called effective permeability, that is the self-averaging property (Boschan and Noetinger, 2012; Noetinger and Gauvier, 1998; Noetinger and Zargar, 2004). It highlights the role of the interaction between the geometry and size of the block and the integral scale that may characterize the internal structure of the aquifer. If these support effects are not accounted for, the simulation of the petrophysical variables can result in distorted and possibly wrong fluid paths. It is therefore crucial to use appropriate geostatistical methods for unstructured grids (and, incidentally, for structured corner point grids in which cell volumes vary significantly, see Bertoncello et al., 2008)).

The simplest approach to account for the support effect is to simulate the properties on a fine-scale grid using a geostatistical model defined at measurement-support. In a second step, this grid is overlaid by the unstructured grid of interest. The values on the large cells are computed by averaging the values contained in this cell on the fine scale grid (e.g., Caumon et al., 2005; Durlofsky, 2005). However, modeling the domain of interest at measurement-support (typically a few cubic centimeter for petrophysical core measurements) is particularly challenging. To bypass this problem, a solution is to simulate directly the property on the unstructured grid accounting for support effect. But a second challenge arises because permeability is not additive. A property is said to be additive (porosity for example) when the averaging can be carried out using a (weighted) arithmetic mean. This is possible because quantities of the underlying physical property can be added together in a meaningful manner. For example, volumes of pores in two parts of a sample can be added together to produce the total volume of pores in the sample.

For additive properties, some direct geostatistical simulation methods accounting for support effect exist. They are based on the analytical or numerical estimation of the covariance between data points and cells of different support and geometries. Among the existing methods, the Discrete Gaussian Model or DGM (Matheron, 1976) has been investigated and implemented by Emery and Ortiz (2011) and Zaytsev et al. (2015). Another possible method is the Direct Sequential Simulation or DSS (Deutsch et al., 2002). In all cases, the covariance values vary depending on the cell's volumes. However, permeability is not additive: one cannot simply average it using simple summations. Instead, a pressure solver or some approximations must be used (Renard and De Marsily, 1997; Manchuk et al., 2012; Khan and Dawson, 2004) and the geostatistical methods cited above cannot be directly applied.

The aim of this paper is, therefore, to introduce a new approach allowing to simulate permeability directly on any unstructured grid, accounting for the support effect and avoiding the use of an underlying fine grid. The method transforms the permeability into an additive property using power averaging with local exponents which depend on the geometry and size of the cells. The exponents are estimated using a lim-

ited set of numerical experiments and an experimental design approach. The simulation process then relies on a Spectral Turning Bands approach (Mantoglou and Wilson, 1982; Emery et al., 2016). This overall strategy permits to generate efficiently a set of realizations. In Section 3, the applicability of the method is illustrated on a synthetic case including tracer test simulations.

## 2. Methodology

The general workflow of the proposed approach is illustrated in Fig. 1. The details are provided in the following sections. Let us first introduce the overall approach. The first step is the computation of the size and aspect ratio of each cell of the unstructured grid (Fig. 1a.). The second step consists of a series of numerical experiments allowing to estimate, for any cell geometry, the exponent $\omega$ that would minimize the error between the equivalent conductivity estimated using a power average with the $\omega$ exponent and the equivalent conductivity estimated numerically. Since it is not possible to run numerical experiments for all the cell sizes and geometries, we select some representative cell dimensions and construct a response surface (Fig. 1b.). This is done once at the beginning of the procedure. Then, several points are randomly placed inside each cell and, for every realization, the permeability is simulated on these points using the Spectral Turning Bands (STB) approach (Fig. 1c.). The power averaging technique is then applied using the value of $\omega$ derived from the response surface to obtain the equivalent permeability of the cell (Fig. 1d.).

The main strength of this workflow is that the averaging exponents $\omega$ are estimated once for a given problem constrained by a permeability distribution, a variogram and an unstructured grid. After that, several realizations of permeability fields can be obtained efficiently using spectral turning bands.

### 2.1. Power averaging

The proposed method is based on power averaging (Deutsch et al., 2002; Journel et al., 1986; Ababou, 1996; Desbarats, 1992) to estimate the equivalent permeability $K^{eq}$ of a cell:

$$K^{eq} = \left( \frac{1}{n} \sum_{i=1}^{n} k_i^{\omega} \right)^{\frac{1}{\omega}}, \quad (1)$$

with $n$ representing the number of points in a block, $k_i$ the permeability value at point $i$ and $\omega$ the averaging exponent. The exponent $\omega$ belongs to the interval $[-1, 1]$, with $\omega = 1$ representing an arithmetic mean, $\omega = -1$ a harmonic mean and $\lim_{\omega \to 0} \dots$ is the geometric mean given by $\exp(\frac{1}{n} \sum_{i=1}^{n} \log(k_i))$. The interest of this approach is that the permeability raised at the exponent $\omega$ becomes an additive variable. However, the main question is then to estimate the value of $\omega$ for a certain configuration of the fine scale permeability. Indeed, it is expected that it may depend on the type of spatial distributions and on the size and geometry of the cells.

Therefore, let us first summarize some theoretical and experimental results from the literature. Several studies have given estimates of the value of $\omega$ for some specific distributions of permeability. In two dimensions, when the permeability $k$ and its inverse ($h = 1/k$) have the same probability distributions, Matheron (1966, 1967, 1968) proves that the equivalent permeability is the geometric mean if the statistics of the spatial distributions of $k$ and $h$ are invariant by rotation, i.e. $\omega = 0$. A case that satisfies the previously mentioned conditions is a 2D isotropic medium with a log normal permeability distribution. In three dimensions, Noetinger (1994) proposes an approximation for this specific case: $\omega = 1/3$. The validity of this formula up to the fourth order (small perturbation) has been demonstrated by Dagan (1993). However, De Wit (1995) and Abramovich and Indelman (1995) later showed using a sixth order development that in three dimensions it is not strictly valid and would necessitate a correction.
<figure><img src="image_4.png" /><figcaption>Fig. 1. General workflow for direct simulation of permeability on unstructured grids. For each unstructured cell, a) we find equivalent dimensions. b) We use these dimensions to constrain a response surface of ω calibrated through upscaling experiments. Then, c) we fill the cells with several integration points and simulated permeability on them. Finally, d) we apply power averaging using the ω estimated in step b).</figcaption></figure>

For anisotropic media, Ababou (1996); Desbarats (1992); Duquerroix et al. (1993); Kruel Romeu (1994); Noetinger and Haas (1996) suggest a simple analytical form. Let $\lambda = L_V/L_H$ be the geometrical anisotropy depending on the ratio of the variogram ranges and $\kappa = k_V/k_H$ be the permeability anisotropy ratio. Changing variables in the Laplace Darcy equation allows to define a global anisotropy ratio through $\alpha = \sqrt{\kappa/\lambda}$. The proposed value of ω is:

$$ \omega(\alpha) = \frac{\arctan(\alpha)}{\pi - \arctan(\alpha)} \qquad (2) $$

Massonnat (2009) introduces $\omega_H$ for the horizontal permeability and $\omega_V$ for the vertical permeability:

$$ \omega_H(\alpha) = \frac{\arctan(\alpha)}{\pi - \arctan(\alpha)} \qquad (3) $$

and

$$ \omega_V = -2 \omega_H + 1 \qquad (4) $$

For isotropic random media, horizontal permeability is approximately equal to vertical permeability with $\omega = 1/3$ as discussed above. Furthermore, to account for the finite size of the cells and non-ergodic effects (the integral scale is not necessarily much smaller than the cell size), Massonnat (2009) introduced two additional coefficients $\epsilon_H$ and $\epsilon_V$ to estimate the global anisotropy ratio:

$$ \alpha = \frac{L_H}{L_V} \sqrt{\frac{k_V}{k_H}} \epsilon_H \epsilon_V \qquad (5) $$

The parameter $\kappa = k_V/k_H$ can be estimated from field measurements, but this is rarely done. A table of values has been proposed for turbiditic facies (Wigniolle and Massonnat, 2013). However, the main difficulty revolves around the dependency of ω to the coefficients $\epsilon_H$ and $\epsilon_V$, that are hardly possible to estimate without empirical curves. Godoy et al. (2018) showed numerically how the exponent is influenced by the block size. The power average was recently revisited by Liao et al. (2020) and references therein. The conclusion from this brief literature review is that numerical experiments must be conducted to estimate ω, as there is no simple analytical solution. In the following section, we will present how the values of ω are defined for all unstructured grid cells, but before explaining that aspect we introduce how we compute an equivalent geometry for the cells of the unstructured grids.

## 2.2. Characterizing the unstructured cells

A preliminary step is to estimate the geometries of the unstructured cells. Indeed, ω depends on the cell's sizes. Two separate sets of axes
<figure><img src="image_5.png" /><figcaption>Fig. 2. Nomenclature of the two sets of axes used in the method. XYZ axes are defined by the grid bounding box orientation. The Mmt axes follow the anisotropy axes of the variogram. The unstructured grid will be considered in Mmt space through projections of the cells on the axes.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Definition of the horizontal aspect ratio H. The distance ΔM_u is obtained by projecting the cell's vertices on the M axis and taking the length between the two extremal points. The same process is applied for the m axis.</figcaption></figure>

are defined (Fig. 2). The XYZ space is characterized by the orientation of the unstructured grid, i.e. the orientation of the bounding box of the grid. The Mmt space is defined following the axis M corresponding to the maximum range of the variogram model in the stratigraphic plane, the axis m corresponding to the minimum range of the variogram model in the stratigraphic plane (perpendicular to M) and the vertical axis t. In this work, grids are 2.5D, which brings Z ≡ t. It is in this depositional space, Mmt, that unstructured cell sizes are defined.

To define the cell sizes, a horizontal aspect ratio H along the Mm axes is defined. It is computed for every cell: all vertices are projected along M and m (Fig. 3). Cell's lengths along axes, ΔM_u and Δm_u, are defined by the distance between the furthest projected vertices on the considered axis. The aspect ratio H is obtained by dividing the projected length along M by the one along m: H = ΔM_u/Δm_u.

For each unstructured cell, an equivalent regular hexahedral cell is defined through a simple system of equations, so that the volume of the

regular cell V_r is equal to the volume of the unstructured cell V_u and that both cells have the same aspect ratio H and thicknesses Δt_r and Δt_u, respectively:

$$ 
\left\{ 
\begin{array}{l} 
V_r = \Delta M_r \times \Delta m_r \times \Delta t_r = V_u \\ 
H = \Delta M_u / \Delta m_r = \Delta M_r / \Delta m_r \\ 
\Delta t_u = \Delta t_r 
\end{array} 
\right. 
\qquad (6) 
$$

ΔM_u, Δm_u, Δt_u and V_u being known, the dimensions of the equivalent rectangular cell are obtained for every cell of the unstructured grid:

$$ 
\left\{ 
\begin{array}{l} 
\Delta M_r = \sqrt{\frac{H \times V_u}{\Delta t_u}} \\ 
\Delta m_r = \sqrt{\frac{V_u}{H \times \Delta t_u}} \\ 
\Delta t_r = \Delta t_u 
\end{array} 
\right. 
\qquad (7) 
$$

For a non-extruded unstructured grid, the third equation of the system should be replaced by a vertical aspect ratio V = (ΔM, Δm)/Δt. Describing the unstructured cells in Mmt space means that the permeability tensor will be considered in this space. For the flow simulators, however, permeability has to be given in XYZ space. An additional step at the end of the workflow has to be implemented to rotate the permeability matrix from Mmt space to XYZ. For simplicity, however, the remainder of this paper considers that Mmt axes are aligned to XYZ axes, coherently with the chosen horizontally isotropic grids and variograms used for the simulations.

### 2.3. Generating the response surface of omega

In the method presented hereafter, ω is estimated for each cell using a response surface methodology. The explanatory variables of the response surface are the different parameters having an influence on the value of ω. They characterize each cell of the unstructured grid.

Following Noetinger and Haas (1996) and Massonnat (2009), experiments show that three main parameters influence ω: variogram ranges, proportions of facies, and cell geometry. In this paper, the focus is put on cell sizes and geometry and we neglect for the moment the presence of multiple facies to keep the response surfaces simple and tractable.

To obtain the response surface, the first step consists in choosing N points in the parameter space where the response has to be evaluated. The dimension D of the parameter space depends on the type of unstructured grid (Fig. 4a.):

* D = 2 for a 2D unstructured grid: taking into account cell's sizes along M and m axes.

* D = 2 for a 3D unstructured grid with constant cell's size on t axis: taking into account cell's sizes along M and m.

* D = 2 for a 3D horizontally isotropic grid: taking into account the cell's sizes along t axis and the average of the cell's sizes along M and m axes.

* D = 3 for a 3D unstructured grid: taking into account cell's sizes along M, m and t axes.

For each of the D dimensions, the minimum and maximum of the sizes of the equivalent rectangular cells (ΔM_r, Δm_r and Δt_r) is computed to limit the ranges of possible cell sizes. A set of P points is then randomly and uniformly placed in the parameter space. Each point corresponds to a D-dimensional vector of the form (ΔM'_r, Δm'_r, Δt'_r). The number of points P should be large enough to cover the entire space (Fig. 4b.). Among those P points, N are chosen using the Wootton, Sergent, Phan-Tan-Luu (WSP) (Sergent, 1989; Sergent et al., 1997a; 1997b) space-filling design technique: from a set of candidate points, well-distributed experiments are selected following an iterative process depending on an exclusion sphere (Santiago et al., 2012)(Fig. 4c.).

A value of ω is then estimated for each of the N points of the parameter space: for each point having D coordinates representing cell sizes, e.g., ΔM_r^{exp}, Δm_r^{exp} and Δt_r^{exp} if D = 3, we use a fine regular grid aligned
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** a) Example of parameter space for $D = 2$ and $D = 3$. b) Placing *P* random points to cover the entire space. c) Choosing *N* experiments uniformly placed in the space.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Definition of the structured grid on which we perform the experiments. We take into account $\Delta M_r^{smallest}$, $\Delta m_r^{smallest}$, $\Delta t_r^{smallest}$ and $\Delta M_r^{biggest}$, $\Delta m_r^{biggest}$, $\Delta t_r^{biggest}$, which are the dimensions of the smallest and biggest unstructured cells.</figcaption></figure>

on the *Mmt* axes to simulate several realizations of fine-scale permeability. These fine values are then upscaled with a pressure solver using upscaling ratio corresponding to the coordinates of the point treated, i.e. we create upscaled blocks of the size of the equivalent hexahedron of interest. It gives a distribution of reference upscaled permeability values and we find the $\omega$ for which the power averaging distribution best fits the reference. The details of this calibration are given below.

The fine grid for the fine-scale simulations is defined once, using $\Delta M_r^{smallest}$, $\Delta m_r^{smallest}$ and $\Delta t_r^{smallest}$ (respectively $\Delta M_r^{biggest}$, $\Delta m_r^{biggest}$ and $\Delta t_r^{biggest}$) which are the dimensions of the equivalent rectangular cell corresponding to the smallest (respectively the biggest) unstructured cell in the grid (Fig. 5 top). The idea is to create a fine grid covering at least eight equivalent rectangular cells (2 along each axis): the bounding box of this structured grid must have dimensions $G_M^s$, $G_m^s$ and $G_t^s$ such that $G_M^s \ge 2\Delta M_r^{biggest}$, $G_m^s \ge 2\Delta m_r^{biggest}$ and $G_t^s \ge 2\Delta t_r^{biggest}$. The voxels of this grid have dimensions $c_M^s$, $c_m^s$ and $c_t^s$ such that $c_M^s \le \Delta M_r^{smallest}$, $c_m^s \le \Delta m_r^{smallest}$ and $c_t^s \le \Delta t_r^{smallest}$ (Fig. 5 bottom).

Then, each experiment is performed as follows:

* A fine-scale permeability *k* is simulated on the local grid using Spectral Turning Bands (Mantoglou and Wilson, 1982; Emery et al.,

2016), with a distribution and variogram corresponding to the fine scale (Fig. 6a.). This step is repeated *R* times with a different seed for random number generation. The permeability distribution can be of any type including lognormal and beta. This distribution may characterize the values of permeability or the values of log-permeability. In the latter case, the fine-scale permeability values are normalized using a power of 10. The rest of the experiment is then performed on this normalized property.

* Upscaling ratios in each direction are set using cell's sizes corresponding to the current experiment: $U_M = \Delta M_r^{exp}/c_M^s$, $U_m = \Delta m_r^{exp}/c_m^s$ and $U_t = \Delta t_r^{exp}/c_t^s$.

* *R* upscaled permeabilities, called reference permeabilities $K_M^{ref}$, $K_m^{ref}$ and $K_t^{ref}$, are computed using a pressure solver upscaler with $U_M$, $U_m$ and $U_t$ as upscaling ratios (Fig. 6b.). It gives a distribution of $K^{ref}$ values in each direction. The upscaler (Jaquet et al., 2000) uses a local scheme and solves the steady state flow equation using a finite volume technique. It uses a Two Point Flux given by Approximation (TPFA) scheme. Two options are available: the type of boundary conditions around the domain (it can either be of permeameter-type or linearly varying heads type, see definition and discussion in Renard and De Marsily (1997)); and the formula used to compute
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** a) Simulation of permeability on the fine grid. This step is performed *R* times. b) Numerical upscaling using the experiment cell sizes to calculate the upscaling ratio. This step is also performed *R* times. c) Optimizing the global error to the reference distribution to obtain an *ω* value.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Examples of *ω* calculation: optimization of the error through a parabola fitting or golden number algorithm.</figcaption></figure>

the transmissibility between two cells (it can either be an arithmetic, geometric or harmonic mean). The resulting linear system is solved using a conjugate gradient solver.

* The values of $\omega_M$, $\omega_m$ and $\omega_t$ for this experiment are determined by minimizing the global error between $K^{\text{ref}}$ and the power average $K^{\text{eq}}$, i.e. minimizing the sum of the error for the $n$ upscaled blocks (Fig. 6c.):

$$ \sum_{\text{block}=1}^{n} \| \log(K_{\text{block}}^{\text{ref}}) - \log(K_{\text{block}}^{\text{eq}}) \| \quad (8) $$

In practice, determining $\omega$ is done using an optimization method that first tries to fit a parabola equation to the error function and find its minimum, and, if not successful, uses a golden search method. The error function can have various shapes, going from a parabola to a simple monotonous curve (Fig. 7).

Typically, the number *R* of fine-scale permeability realizations must be large enough so that the reference distribution $K^{\text{ref}}$ has a statistically sufficient number of values to identify a robust value for $\omega$. In the case of the biggest possible equivalent hexahedron (of size $\Delta M^{\text{biggest}}$, $\Delta m_r^{\text{biggest}}$ and $\Delta t_r^{\text{biggest}}$), we obtain eight upscaled values for one realization: the fine grid has been defined such that eight biggest hexahedra can be overlaid on it. It means that in this case, $\omega$ will be fitted on $n = 8R$ values, which is the limit to take into account when choosing the number *R*. In parallel, in the case of a small equivalent hexahedron, we will obtain a large number of reference values for each realization, because the fine grid is not adapted to the local problem. This large

number will also be multiplied by *R*, which gives too many unnecessary values to estimate $\omega$ and brings unnecessary computational cost. To fix this problem, if the cell size for the treated experiment is smaller than $\Delta M_r^{\text{biggest}} \times \Delta m_r^{\text{biggest}} \times \Delta t_r^{\text{biggest}}$, only the region of the fine grid of size $2\Delta M_r^{\text{exp}}$, $2\Delta m_r^{\text{exp}}$ and $2\Delta t_r^{\text{exp}}$ is simulated and consequently upscaled. By doing so, $K^{\text{ref}}$ distribution always has $n = 8R$ values to fit the $\omega$ value for each point in parameter space.

Once all numerical experiments have been performed, an interpolation between the points of the parameter space is done using kriging (Fig. 8a. and b.). The correlation function for kriging is assumed stationary and written as follows:

$$ R(x, y) = \exp\left(-\sum_{j=1}^{D} \left(\frac{|x_j - y_j|}{\theta_j}\right)^p\right) \quad (9) $$

with $1 < p \le 2$ and *D* the number of dimensions of the response surface. The hyper-parameters $\theta_j$, characterizing kriging anisotropy, are estimated using the maximum likelihood method. However, to avoid artifacts on surfaces of response in our work, all $\theta_j$ have been set to one in the normalized space. The resulting response surface allows determining the exponent $\omega$ for each cell of the unstructured grid (Fig. 8c.).

As explained previously, in each experiment the upscaling is done in two or three directions resulting in different $K_{\text{ref}}$ (*M*, *m* or *t*). Then, two or three values for $\omega$ are obtained per experiment, one for each direction *D*. In the end, not one but two or three response surfaces of $\omega$ are interpolated and the same amount of omega properties are painted on the unstructured grid.
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** a) All experiment results, i.e. $\omega$, are put in the parameter space. b) A simple kriging allows the generation of a surface of response of $\omega$. c) By projecting cell sizes on the response surface, one $\omega$ is found for every cell of the unstructured grid.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** a) One by one, each unstructured cell is filled with points, the number of points depending on the volume of the cell; b) Using the Spectral Turning Bands, permeability is interpolated on these points; c) Using the $\omega$ of the target cell, the permeability of the cell is obtained using power averaging on the points (followed by a back transform from $K^\omega$ to $K$).</figcaption></figure>

## 2.4. Direct simulation of permeability on the unstructured grid

The next step is to simulate rapidly values on integration points using Spectral Turning Bands (STB) and average them using the estimated $\omega$. The integration points are obtained with a fast generation of quasi-random Sobol' sequences using the method of Antonov and Saleev (1979).

Turning bands were introduced in Matheron (1973) and further developed by many authors including Mantoglou and Wilson (1982), Tompson et al. (1989), Lantuéjoul (2002) and Emery and Lantuéjoul (2006). The general principle of this technique is to reduce a D dimensions simulation problem to a one-dimensional simulation. The first step is to generate L lines $\theta_i$ such that their orientation is uniformly distributed over the unit sphere. A 1D simulation is made along each line using the one-dimensional field covariance $C_\theta$ resulting from the D-dimensional covariance C. Numerous methods for simulating a one-dimensional random field knowing its covariance $C_\theta$ are found in the literature, including spectral approaches used in this work. They can be classified into continuous spectral simulation using cosine functions (Mantoglou and Wilson, 1982; Shinozuka and Jan, 1972); discrete spectral simulation using Fast Fourier Transform (Tompson et al., 1989); and circulant embedding (Dietrich, 1995). From these simulated lines, the simulated value at any point can be obtained by projection and summation. An extension of the Spectral Turning Bands method to non-stationary models is available in Emery and Arroyo (2017).

Once the lines have been simulated and stored, the simulation of each unstructured cell can be made independently from the others. For each cell, a number of random locations to select is calculated depending on the cell's volume (Fig. 9a.). The fine permeability values at these locations are simulated as explained above (Fig. 9b.). The global cell value is obtained by averaging the point values using the power law and the $\omega$ estimated for that cell (Fig. 9c.). Once again, not one but two or three permeability values are obtained, one for each direction M, m and t.

## 3. Applications

The methodology has been tested and showed to provide correct results on simple cases with known analytical solutions (perfectly layered case, isotropic log-normal distribution, etc.). We do not discuss these results here for the sake of brevity.

Instead, we will illustrate the method with a synthetic but realistic example (Fig. 10). We consider a confined aquifer of 2km by 2km with negligible ambient flow. The external boundaries are assumed to be impermeable. The thickness of the aquifer is varying linearly from 40 m in the eastern side to 4 m in the western side. Ten wells are positioned across the aquifer. Two of them are injecting water and seven are pumping (Fig. 15). A tracer is injected in the last well and recovered in the pumping wells. This setup was designed to test the proposed geostatistical simulation technique on a realistic unstructured grid and to compare the results of the tracer simulation with a fine scale simulation (standard approach).

In the following sections, we describe first the simulation of the permeability and then the tracer simulations.

### 3.1. Permeability simulation

The 3D unstructured grid for this case is composed of 18,250 Voronoï cells with a volume between $25 \text{ m}^3$ and $50.10^3 \text{ m}^3$ (Fig. 10a.). The grid is refined around all ten wells. The variogram for both porosity and permeability has two Gaussian structures of equal contribution with ranges of $25 \times 25 \times 2.5 \text{ m}$ and $250 \times 250 \times 25 \text{ m}$ for cell sizes varying between 3 and 90 m horizontally and 0.8 to 8 m vertically (Fig. 10b.). The porosity has a beta distribution of mean 0.15 and standard deviation 0.05 (Fig. 10c.) and the permeability distribution is lognormal with a mean of 100 mD and a standard deviation of 300 mD (Fig. 10d.). The boundary conditions for the upscaler in the experiments are of permeameter-type, i.e. no-flow conditions on the sides.

Due to the horizontal isotropy of the aspect ratio of the unstructured cells, the Mmt axes presented in Section 2.2 are taken aligned to XYZ axes. For the same reason, X axis and Y axis being equivalent, only two response surfaces are presented hereafter: one for the horizontal exponent $\omega_H = \omega_X = \omega_Y$ and one for the vertical exponent $\omega_V$ (Fig. 11). The two axes of the response surfaces represent the horizontal cell sizes $\Delta XY$ (taken as the mean of the cell sizes in X and Y), and the vertical cell sizes $\Delta Z$. A smooth variation of $\omega$ is observed, from 0 to 1 for $\omega_H$ and from -1 to 1 for $\omega_V$.

The log of permeability is assumed to be correlated to point-support porosity with a correlation coefficient of 0.8, which is possible using Spectral Turning Bands. Porosity on the unstructured grid is obtained using the arithmetic mean (Fig. 12) and two permeability (H and V) are simulated using the power averaging formula and $\omega$ properties gener-
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** a) grid for the synthetic case with cell volumes varying with a factor 2.10³, b) double structure variogram used for the simulation of both petrophysical variables. c) beta distribution of porosity, with m=0.15 and σ=0.05. d) lognormal distribution of permeability, with m=100mD and σ=300mD.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Surfaces of response of horizontal and vertical omega and the corresponding properties on the unstructured grid.</figcaption></figure>

<caption>**Fig. 12.** Porosity field on the toy example. Log-permeability will be correlated to this field with a linear coefficient equal to 0.8.</caption>

ated previously (Fig. 13a.). Although both $K_H$ and $K_V$ properties look similar, the cross-plot on Fig. 13b. shows that $K_V$ is globally inferior to $K_H$. The similarity between these results will be discussed in the conclusion.

## 3.2. Tracer tests

Tracer simulations are performed accounting only for advection and neglecting diffusion and physical dispersion. The tracer response computed using the unstructured grid presented above is compared to the response computed using a fine structured grid having the same extension. The aim is to test if the permeability field obtained with the proposed method gives a coherent tracer response. The fine grid has 600 × 600 × 80 cells of horizontal size 3.3 × 3.3 m and vertical size from 0.05m to 0.5m. Using spectral turning bands (STB), it is possible to obtain the same random field on the points placed in the structured and unstructured cells. For the structured one, a single point is taken per cell, while the number of points in each unstructured cells depends on its volume. Then, for unstructured cells, the cell average value is obtained using arithmetic average for porosity and power averaging for permeability (using the $\omega$ value given by the response surface). Twenty realizations of porosity and correlated
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** a) Vertical and horizontal permeability fields on the toy example. b) Cross-plot $K_H$ vs. $K_V$, $K_V$ being globally inferior to $K_H$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14.** a) Simulation of porosity and permeability on the fine structured grid. b) Simulation of porosity and permeability on the coarser unstructured grid. Details are maintained in the most refined zones.</figcaption></figure>

permeability are generated on both grids, one example is displayed on Fig. 14.

For each realization, a passive tracer is injected in a water saturated reservoir having one well injecting the tracer, two wells injecting water and seven production wells (Fig. 15).

The water production rate per well is 120 m³/day and the water (and tracer) injection rate per well is 280 m³/day. We observe the concentration of tracer at the producers for the twenty realizations (Fig. 16).

For individual realizations, some local deviations are observed between the tracer curves of the fine structured reservoir (Fig. 16a.) and the ones of the unstructured grid (Fig. 16b.). However, in terms of uncertainty assessment, the set of tracer curves shows a very consistent behavior for both grids. The curves for structured reservoir are slightly more scattered, which is coherent with the better precision offered by this grid. The Q10, Q50 and Q90 quantile curves (Fig. 16c.) are almost identical, with some minor differences on Q10. A 3D visualization of the tracer's concentration evolution shows similar results on structured and

unstructured grids (Fig. 17). This similarity is further demonstrated by the visualization of the mean of the twenty realizations of the tracer's concentration evolution at time thirteen years, where concentration tendencies are well respected (Fig. 18). As expected for a coarser grid, these results show that the unstructured grid does not preserve the details of the propagation visible on the structured grid, horizontally and vertically. This highlights the importance of the choice of unstructured grid: the errors are localized in areas between wells where the unstructured grid is too coarse. Overall, this comparison demonstrates the consistency of the proposed method as compared to the results obtained with a structured grid. However, the computational times are very different. The tracer test simulations on the structured grid takes 112 minutes in average, while they only take 4 minutes on the unstructured grid.

### 3.3. Computation times

The use of an unstructured grid, while accounting for the support effect, usually requires to generate the geostatistical simulations on a
<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Placement of the wells in the unstructured reservoir (painted with a realization of porosity upscaled arithmetically from points simulated with STB). There are seven producing wells, one injecting tracer and two injecting water. The placement of wells is identical on the fine structured reservoir model.</figcaption></figure>

<caption>**Table 1**</caption>

Computation times for two unstructured grids of different sizes. $T_0$ is the time to generate the response surfaces and $\omega$ properties on the grids, $T_3$ is the time to simulate permeability on points using STB and perform power averaging.

<table><thead><tr><th></th><th>Small grid</th><th>Large grid</th></tr></thead><tbody><tr><td>nb cells</td><td>18,250</td><td>2,031,995</td></tr><tr><td>$T_0$</td><td>22.4 min</td><td>35.4 min</td></tr><tr><td>$T_3$</td><td>11.3 s</td><td>7.5 min</td></tr><tr><td>Total (N=1)</td><td>22.4 min</td><td>43.0 min</td></tr><tr><td>Total (N=10)</td><td>24.3 min</td><td>111.1 min</td></tr></tbody></table>

fine structured grid, and to upscale numerically the permeability (e.g. Aavatsmark et al., 1998; Prévost et al., 1996). For every realization, one needs to redo the upscaling. The methodology presented in this paper is different since it allows direct simulation of the permeability on the un-structured grid, but it shows improved efficiency for multi-realizations.

In particular, the computation time for the proposed method is not proportional to the number of cells in the unstructured grid. Indeed, the response surface computation is generally the most computationally demanding step, so the main factor influencing the time is the range of cell sizes covered. A grid presenting a larger range of cell sizes will require a broader response surface and therefore more experiments. In addition, the fine structured grid used to perform the experiments is defined from the minimum and maximum cell sizes. A higher difference between the two will require a larger fine grid, and hence longer upscaling times.

To evaluate the computation time, different tests have been made. The procedure has been applied on two grids with different characteristics: a small example with 18,250 cells varying from 3 to 90 m horizontally and from 0.8 to 8 m vertically, and a larger example with 2,031,995 cells varying from 6 to 2,653 m horizontally and from 3 to 60 m vertically.

In the following, $T_0$ represents the time required to generate the response surfaces, $T_3 = T_1 + T_2$ represents the time required to simulate the permeability values on points using STB ($T_1$) and perform power averaging ($T_2$). The total simulation time $T$ is:

$$T = T_0 + N T_3 \qquad (10)$$

with $N$ the number of realizations.

For a response surface with fifty experiments, using twenty Intel cores on a Linux machine, the computation times for our Java implementation are provided in Table 1. This table shows that the time required

<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Concentration of tracer at producers curves for twenty realizations on a) the fine structured reservoir and b) the unstructured reservoir. The colours of the curves for a) and b) represent the same realizations. c) Quantiles of the twenty curves for structured (in blue) and unstructured (in red) reservoirs. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

to perform one realization on the large grid is only doubled compared to the time required for the small grid, while the number of cells has been multiplied by a factor 100, showing the efficiency of the method for large grids.

The computation time is also compared with the classical method (fine grid simulation followed by a pressure solver upscaling on a regular grid). The experiment is conducted as follows. The fine grid is the one presented in Section 3.2, with 600x600x80 cells, the coarse grid has 60x60x5 cells (the upscaling ratios are 10x10x16). This coarse grid has 18,000 regular cells, which is close to the 18,250 cells of the unstructured grid.

For one realization, the proposed method takes three minutes more than the classical one (Table 2). However, as soon as one needs to perform at least two realizations, the proposed methodology is faster. Moreover, the response surfaces being saved, it is almost immediate to simulate permeability fields when parameters have not changed, i.e. when the distribution, variogram and grid are the same.
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Visualization of the concentration of tracer in the structured and unstructured reservoirs at a) seven years of simulation and b) thirty years of simulation.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Visualization of a) the mean and b) the standard deviation of the twenty realizations of tracer's concentration evolution at time thirteen years in the structured and unstructured reservoirs.</figcaption></figure>

<caption>**Table 2**</caption>

Computation times for the classical method (fine grid simulation and pressure solver upscaling) and the proposed method. $T_0$ = time to generate the response surfaces, $T_1$ is the time for STB simulation on points and $T_2$ is the time for either power averaging on the unstructured grid or pressure solver upscaling for the structured grid.

<table><thead><tr><td></td><td>Proposed method</td><td>Classical method</td></tr></thead><tbody><tr><td>$T_0$</td><td>22.4 min</td><td>0</td></tr><tr><td>$T_1$</td><td>x</td><td>4.6 min</td></tr><tr><td>$T_2$</td><td>x</td><td>15.1 min</td></tr><tr><td>$T_3 = T_1 + T_2$</td><td>11.3s</td><td>19.8 min</td></tr><tr><td>Total (N=1)</td><td>22.4 min</td><td>19.8 min</td></tr><tr><td>Total (N=2)</td><td>22.8 min</td><td>39.5 min</td></tr><tr><td>Total (N=10)</td><td>24.3 min</td><td>394.7 min</td></tr></tbody></table>

## 4. Discussion and conclusion

This paper proposes a new workflow allowing to simulate directly permeability fields on unstructured grids accounting for support effects.

The main originality of the method is that it allows bypassing the use of a potentially memory and time consuming fine grid and the repeated use of local upscaling on every fine-scale realization. On horizontally isotropic simple cases, the applications are encouraging, with coarse flow simulations results close to a fine grid taken as reference. The proposed method can be used as an extension to previous methods that allowed to generate geostatistical simulation of additive variables directly on an unstructured grid (e.g. Zaytsev et al., 2015; Deutsch et al., 2002; Emery and Arroyo, 2017).
The methodology raises some questions that are detailed below. The obtained results provide some clues to answer these questions and open interesting discussion items.

The first one relates to the numerical upscaling method used in the experiments as a reference. Numerical upscaling is more general than analytical solutions, but it is strongly dependent on a few parameters such as the type of boundary conditions, or the choice of the averaging technique used for transmissibility computation. Three main boundary conditions are used in general: permeameter-type, linearly-varying head, and periodic. In this paper, we presented results that were obtained with permeameter-type conditions. However, the same cases could be studied with linearly-varying head conditions. To select which method is the most adequate, additional research should be conducted. The general methodology would, however, not be affected.

The second one is whether the power averaging formula is sufficient to capture the details of the spatial complexity of the permeability fields. Previous numerical experiments (e.g. Renard et al., 2000) have shown for example a broad dispersion of the equivalent permeabilities around power averages. However, and rather surprisingly, the numerical experiments conducted in this paper show that the method is robust for the studied class of random fields. In addition, one can note that during the optimization step when the value of $\omega$ is estimated from a set of reference upscaled permeabilities, we store the remaining error. This means that, even if we did not encounter this problem yet, the methodology includes a step to identify situations in which the method could have difficulties. In these situations, one could apply a quantile-to-quantile correction of permeability, replacing the permeability in cells flagged as “bad” by the corresponding reference, numerically upscaled, permeability.

The global consistency of the approach can be checked by using the same simulator at the finest and coarse scales. It is well known from Romeu and Noetinger (1995) that any discretized model can yield strongly biased results as soon as the grid block size is on the same order as the underlying correlation scale. This bias is due to the underlying numerical scheme that weights the conductivity between grid-blocks by means of some averaging formula that may lead to underestimated values for the overall effective conductivity. That is observed in the popular case of the harmonic averaging formula employed by most popular commercial simulators. In the common practice, the mesh of the working model is built by engineers in order to get the best compromise: ensuring a global accuracy at the lowest numerical cost. Close to the wells, very fine grid blocks are employed, while coarse blocks are used far from the wells. The bias issue may thus occur in the transition zone with grid block sizes that may be compared to the size of the heterogeneities. Such a zone may be expected to be of a rather small size compared to the overall characteristic sizes of the model (distances between wells, size of the reservoir). Similar issues were addressed by Preux (2016) in order to check the location at which upscaling should be carried out. A posteriori estimators may be used to postprocess the solutions in order to indicate zones to be refined, see Gratien, Jean-Marc et al. (2016) and references therein. Regarding the boundary condition issue, it can be observed that the effective conductivity of large blocks becomes rather independent on the boundary conditions, see Colecchio et al. (2020) and references therein.

Finally, we want to note that the comparison of the computing time for the standard upscaling approach and the proposed one is probably in favor of the standard one. Indeed, the standard approach used in this paper assumes that the upscaling is done on a regular grid. Computations are easier than if they were done for any possible geometry. However we showed that the proposed method becomes faster when the number of realizations increases, so we consider that the all cost of designing the experiment, running them and interpolating the $\omega$ values is rapidly compensated.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

The authors would like to thank Total for funding this research, and especially the research project Reservoir Earth Modeling for the financial and technical support through the help of Jean-Paul Rolando, Gérard Massonnat and David Ledez. We also want to thank the Association Nationale de la Recherche et de la Technologie (ANRT) which funded this PhD under the convention CIFRE 2017/0222.

## References

Aavatsmark, I., Barkve, T., Bøe, O., Mannseth, T., 1998. Discretization on unstructured grids for inhomogeneous, anisotropic media. Part I: derivation of the methods. SIAM J. Sci. Comput. 19 (5), 1700–1716.

Ababou, R., 1996. Random porous media flow on large 3-d grids: numerics, performance, and application to homogenization. In: Environmental Studies. Springer, pp. 1–25.

Abramovich, B., Indelman, P., 1995. Effective permittivity of log-normal isotropic random media. J. Phys. A 28, 693–700.

Akhurst, M., Hannis, S.D., Quinn, M.F., Shi, J.-Q., Koenen, M., Delprat-Jannaud, F., Lecomte, J.-C., Bossie-Codreanu, D., Nagy, S., Klimkowski, Ł., et al., 2015. Risk assessment-led characterisation of the Sitechar UK North Sea site for the geological storage of CO₂. Oil Gas Sci. Technol. 70 (4), 567–586.

Antonov, I., Saleev, V., 1979. An economic method of computing $lp_1$-sequences. USSR Comput. Math. Math. Phys. 19 (1), 252–256.

Bertoncello, A., Caers, J., Biver, P., Caamon, G., 2008. Geostatistics on stratigraphic grids. In: Ortiz, J., Emery, X. (Eds.), Eighth Geostatistical Geostatistics Congress. Gecamin, Ltd, Santiago, Chile, Chile, pp. 677–686.

Blessent, D., Therrien, R., Gable, C.W., 2011. Large-scale numerical simulation of groundwater flow and solute transport in discretely-fractured crystalline bedrock. Adv. Water Resour. 34 (12), 1539–1552. https://doi.org/10.1016/j.adwatwres.2011.09.008.

Boschan, A., Noetinger, B., 2012. Scale dependence of effective hydraulic conductivity distributions in 3d heterogeneous media: a numerical study. Transp. Porous Media 94 (1), 101–121.

Caumon, G., Grosse, O., Mallet, J.-L., 2005. High Resolution Geostatistics on Coarse Unstructured Flow Grids. In: Leuangthong, O., Deutsch, C.V. (Eds.), Geostatistics Banff 2004. Springer Netherlands, Dordrecht, pp. 703–712.

Chilès, J.-P., Delfiner, P., 2012. Geostatistics: Modeling Spatial Uncertainty, 2nd edition. Wiley.

Colecchio, I., Boschan, A., Otero, A.D., Noetinger, B., 2020. On the multiscale characterization of effective hydraulic conductivity in random heterogeneous media: a historical survey and some new perspectives. Adv. Water Resour. 103594.

Dagan, G., 1993. Higher-order correction of effective permeability of heterogeneous isotropic formations of lognormal conductivity distribution. Transp. Porous Media 12, 279–290.

De Marsily, G., Delay, F., Teles, M., Schafmeister, T., 1998. Some current methods to represent the heterogeneity of natural media in hydrogeology. Hydrogeol. J. 6 (1), 115–130.

De Wit, A., 1995. Correlation structure dependence of the effective permeability of heterogeneous porous media. Phys. Fluids 7 (11), 2553–2562.

Desbarats, A.J., 1992. Spatial averaging of hydraulic conductivity in three-dimensional heterogeneous porous media. Math. Geol. 24 (3), 249–267.

Deutsch, C., Journel, A.G., 1992. GSLIB - Geostatistical Software Library and User's Guide. Oxford University Press.

Deutsch, C., Tran, T.T., Pyrcz, M.J., 2002. Geostatistical assignment of reservoir properties on unstructured grids. In: Proceedings of the SPE Annual Technical conference and exhibition, San Antonio, Texas.

Dietrich, C.R., 1995. A simple and efficient space domain implementation of the turning bands method. Water Resour. Res. 31 (1), 147–156.

Duqueroix, J., Lemouzy, P., Noetinger, B., Romeu, R., et al., 1993. Influence of the permeability anisotropy ratio on large-scale properties of heterogeneous reservoirs. In: SPE Annual Technical Conference and Exhibition.

Durlofsky, L., 2005. Upscaling and gridding of fine scale geological models for flow simulation. In: Proceedings of the 8th International Forum on Reservoir Simulation, Stresa, Italy.

Emery, X., Arroyo, D., 2017. On a continuous spectral algorithm for simulating non-stationary gaussian random fields. Stoch. Environ. Res. Risk. Assess.

Emery, X., Arroyo, D., Porcu, E., 2016. An improved spectral turning-bands algorithm for simulating stationary vector gaussian random fields. Stoch. Environ. Res. Risk Assess. 30 (7), 1863–1873.

Emery, X., Lantuëjoul, C., 2006. Tbsim: a computer program for conditional simulation of three-dimensional gaussian random fields via the turning bands method. Comput. Geosci. 32 (10), 1615–1628.
Emery, X., Ortiz, J., 2011. Two approaches to direct block-support conditional co-simulation. Comput. Geosci. 37 (8), 1015-1025.

Godoy, V.A., Zuquette, L.V., Gómez-Hernández, J.J., 2018. Stochastic analysis of three-dimensional hydraulic conductivity upscaling in a heterogeneous tropical soil. Comput. Geotech. 100, 174-187.

Goovaerts, P., 1997. Geostatistics for Natural Resources Evaluation. Oxford University Press.

Gratien-Jean-Marc, Ricois, Olivier, Yousef, Soleiman, 2016. Reservoir simulator runtime enhancement based on a posteriori error estimation techniques. Oil Gas Sci. Technol. 71 (5), 59. https://doi.org/10.2516/ogst/2016009.

Jaquet, O., Namar, R., Siegel, P., Unpublished reference. UPSCALER: Guide Technique & D'utilisation, Memorandum Upscaler 009/004. Report. In2Earth.

Journel, A.G., Deutsch, C., Desbarats, A.J., 1986. Power averaging for block effective permeability. In: Proceedings of the SPE California Regional Meeting, Oakland, California.

Khan, S.A., Dawson, A.G., 2004. Method of Upscaling Permeability for Unstructured Grids. US Patent 6,826,520 B1

Kruel Romeu, R., 1994. Ecoeulement en Milieux Heterogenes: Prise de Moyenne de Permeabilité en Regimes Permanent et Transitoire. Paris 6 Ph.D. thesis.

Lantuéjoul, C., 2002. Geostatistical Simulation: Models and Algorithms. Springer-Verlag Berlin Heidelberg.

Liao, Q., Lei, G., Wei, Z., Zhang, D., Patil, S., 2020. Efficient analytical upscaling method for elliptic equations in three-dimensional heterogeneous anisotropic media. J. Hydrol. 124560.

Manchuk, J.G., Mlacnik, M.J., Deutsch, C., 2012. Upscaling permeability to unstructured grids using the multipoint flux approximation. Pet. Geosci. 18 (2), 239-248.

Mantoglou, A., Wilson, J.L., 1982. The turning bands method for simulation of random fields using line generation by a spectral method. Water Resour. Res. 18 (5), 1379-1394.

Massonnat, G., 2009. Method, Program and Computer System for Scaling the Modelling Data for a Hydrocarbon Deposit. Patent WO 034253A1

Matheron, G., 1966. Structure et composition des perméabilités. Revue de l'Institut Français du Pétrole XXI, 564-580.

Matheron, G., 1967. Composition des perméabilités en milieu poreux hétérogène, méthode de schwydler et règles de pondération. Revue de l'Institut Français du Pétrole XXII (3), 443-466.

Matheron, G., 1968. Composition des perméabilités en milieu poreux hétérogène: critique de la règle de pondération géométrique. Revue de l'Institut Français du Pétrole XXIII, 201-218.

Matheron, G., 1973. The intrinsic random functions and their applications. Adv. Appl. Probab. 5, 439-468.

Matheron, G., 1976. Forecasting block grade distributions: the transfer functions. In: Guarascio, M., David, M., Huijbregts, C. (Eds.), Advanced Geostatistics in the Mining Industry. Springer Netherlands, pp. 237-251.

Merland, R., Caumon, G., Lévy, B., Collon-Drouaillet, P., 2014. Voronoi grids conforming to 3d structural features. Comput. Geosci. 18, 373-383.

Michael, K., Golab, A., Shulakova, V., Ennis-King, J., Allison, G., Sharma, S., Aiken, T., 2010. Geological storage of CO₂ in saline aquifers: a review of the experience from existing storage operations. Int. J. Greenh. Gas Control 4, 659-667.

Mlacnik, M.J., Harrer, A.W., Heinemann, Z.E., 2003. Locally streamline-pressure-potential-based pebi grids. In: Proceedings of the SPE Reservoir Symposium, Houston, Texas.

Noetinger, B., 1994. The effective permeability of a heterogeneous porous medium. Transp. Porous Media 15, 99-127.

Nøetinger, B., Gautier, Y., 1998. Use of the fourier-laplace transform and of diagrammatical methods to interpret pumping tests in heterogeneous reservoirs. Adv. Water Resour. 21 (7), 581-590.

Noetinger, B., Haas, A., 1996. Permeability averaging for well tests in 3d stochastic reservoir models. In: Proceedings of the SPE Annual Technical Conference and Exhibition, Denver, Colorado.

Noetinger, B., Zargar, G., 2004. Multiscale description and upscaling of fluid flow in subsurface reservoirs. Oil Gas Sci. Technol. 59 (2), 119-139.

Preux, C., 2016. About the use of quality indicators to reduce information loss when performing upscaling. Oil Gas Sci. Technol. 71 (1), 7.

Prévost, M., Lepage, F., Durlofsky, L., Mallet, J.-L., 1996. Unstructured 3d gridding and upscaling for coarse modelling of geometrically complex reservoirs. Pet. Geosci. 11, 339-345.

Quinlivan, P., Batten, A., Wibowo, M., Hinchliffe, S., Rahayu, D., Doria, I., Yahmadi, A., Tondang, H.Y.T., 2015. Assessing geothermal tariffs in the face of uncertainty, a probabilistic approach. In: Proceedings World Geothermal Congress.

Renard, P., De Marsily, G., 1997. Calculating equivalent permeability: a review. Adv. Water Resour. 20 (5-6), 253-278.

Renard, P., Le Loc'h, G., Ledoux, E., de Marsily, G., Mackay, R., 2000. A fast algorithm for the estimation of the equivalent hydraulic conductivity of heterogeneous media. Water Resour. Res. 36 (12), 3567-3580.

Romeu, R., Noetinger, B., 1995. Calculation of internodal transmissivities in finite difference models of flow in heterogeneous porous media. Water Resour. Res. 31 (4), 943-959.

Santiago, J., Claeys-Bruno, M., Sergent, M., 2012. Construction of space-filling designs using wsp algorithm for high dimensional spaces. Chemometr. Intell. Lab. Syst. 113, 26-31.

Sergent, M., 1989. Contribution de la Méthodologie de la Recherche Expérimentale à L'élaboration de Matrices Uniformes: Application Aux Effets de Solvants et de Substituants. Aix-Marseille 3 Ph.D. thesis.

Sergent, M., Phan Tan Luu, R., Elguero, J., 1997. Statistical analysis of solvent scales. Part 1. An. Quim. 93 (2), 71-75.

Sergent, M., Phan Tan Luu, R., Elguero, J., 1997. Statistical analysis of solvent scales. Part 2. An. Quim. 93 (5), 295-300.

Shinozuka, M., Jan, C.M., 1972. Digital simulation of random processes and its applications. J. Sound Vib. 25 (1), 111-128.

Tidwell, V.C., Wilson, J.L., 1997. Laboratory method for investigating permeability upscaling. Water Resour. Res. 33 (7), 1607-1616.

Tompson, A.F.B., Ababou, R., Gehlar, L.W., 1989. Implementation of the three-dimensional turning bands random field generator. Water Resour. Res. 25 (10), 2227-2243.

Vogt, C., Mottaghy, D., Wolf, A., Rath, V., Pechnig, R., Clauser, C., 2010. Reducing temperature uncertainties by stochastic geothermal reservoir modelling. Geophys. J. Int. 181 (1), 321-333.

Wigniolle, E., Massonnat, G., 2013. Use of well test interpretations to constrain the small scale petrophysical anisotropy ($k_v/k_h$). In: Proceedings of the SPE Reservoir Characterisation and Simulation Conference and Exhibition, Abu Dhabi, UAE.

Witter, J.B., Trainor-Guitton, W.J., Siler, D.L., 2019. Uncertainty and risk evaluation during the exploration stage of geothermal development: a review. Geothermics 78, 233-242. https://doi.org/10.1016/j.geothermics.2018.12.011.

Zaytsev, V., Biver, P., Wackernagel, H., Allard, D., 2015. Change-of-support models on irregular grids for geostatistical simulation. Math. Geosci.