Direct simulation of non-additive properties on unstructured grids

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


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 2 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 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 selfaveraging property ( Boschan and Noetinger, 2012;N oe tinger and Gautier, 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, incidently, 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 Ortiz (2011) andZaytsev 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.

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. 1 a.). The second step consists of a series of numerical experiments allowing to estimate, for any cell geometry, the exponent that would minimize the error between the equivalent conductivity estimated using a power average with the 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. 1 b.). 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. 1 c.). The power averaging technique is then applied using the value of derived from the response surface to obtain the equivalent permeability of the cell ( Fig. 1 d.).
The main strength of this workflow is that the averaging exponents 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.

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: with n representing the number of points in a block, k i the permeability value at point i and the averaging exponent. The exponent belongs to the interval [−1 , 1] , with = 1 representing an arithmetic mean, = -1 a harmonic mean and lim →0 ⋯ is the geometric mean given by . The interest of this approach is that the permeability raised at the exponent becomes an additive variable. However, the main question is then to estimate the value of 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 for some specific distributions of permeability. In two dimensions, when the permeability k and its inverse ( ℎ = 1∕ ) have the same probability distributions, Matheron (1966Matheron ( , 1967Matheron ( , 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. = 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: = 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. For anisotropic media, Ababou (1996) ; Desbarats (1992) ; Duquerroix et al. (1993) ; Kruel Romeu (1994) ; Noetinger and Haas (1996) suggest a simple analytical form. Let = ∕ be the geometrical anisotropy depending on the ratio of the variogram ranges and = ∕ be the permeability anisotropy ratio. Changing variables in the Laplace Darcy equation allows to define a global anisotropy ratio through = √ ∕ . The proposed value of is: Massonnat (2009) introduces H for the horizontal permeability and V for the vertical permeability: and For isotropic random media, horizontal permeability is approximately equal to vertical permeability with = 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 H and V to estimate the global anisotropy ratio: The parameter = ∕ 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 H and 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.

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  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 : = Δ ∕Δ . 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: Δ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: For a non-extruded unstructured grid, the third equation of the system should be replaced by a vertical aspect ratio = ⟨Δ , Δ ⟩∕Δ . 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.

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. 4  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 ( Δ , Δ , Δ ). The number of points P should be large enough to cover the entire space ( Fig. 4 b.). Among those P points, N are chosen using the Wootton, Sergent, Phan-Tan-Luu (WSP) ( Sergent, 1989;Sergent et al., 1997a;1997b ) spacefilling 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. 4 c.).
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., Δ exp , Δ exp and Δ exp if = 3 , we use a fine regular grid aligned  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 for which the power averaging distribution best fits the reference. The details of this calibration are given bellow.
The fine grid for the fine-scale simulations is defined once, using Δ , Δ and Δ (respectively Δ , Δ and Δ ) 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 , and such that ≥ 2Δ , ≥ 2Δ and ≥ 2Δ . The voxels of this grid have dimensions , and such that ≤ Δ , ≤ Δ and ≤ Δ ( 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. 6 a.). 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: = Δ exp ∕ , = Δ exp ∕ and = Δ exp ∕ .
• R upscaled permeabilities, called reference permeabilities , and , are computed using a pressure solver upscaler with U M , U m and U t as upscaling ratios ( Fig. 6 b.). It gives a distribution of K ref values in each direction. The upscaler ( Jaquet et al., 0000 ) 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 permeametertype or linearly varying heads type, see definition and discussion in Renard and De Marsily (1997) ); and the formula used to compute 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. 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 M , m and t for this experiment are determined by minimizing the global error between K ref and the power average K eq , i.e. minimizing the sum of the error for the n upscaled blocks ( Fig. 6 c.): In practice, determining 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 ref has a statistically sufficient number of values to identify a robust value for . In the case of the biggest possible equivalent hexahedron (of size Δ , Δ and Δ ), 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, will be fitted on = 8 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 and brings unnecessary computational cost. To fix this problem, if the cell size for the treated experiment is smaller than Δ × Δ × Δ , only the region of the fine grid of size 2Δ exp , 2Δ exp and 2Δ exp is simulated and consequently upscaled. By doing so, K ref distribution always has = 8 values to fit the 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. 8 a. and b.). The correlation function for kriging is assumed stationary and written as follows: with 1 < p ≤ 2 and D the number of dimensions of the response surface. The hyper-parameters j , characterizing kriging anisotropy, are estimated using the maximum likelihood method. However, to avoid artifacts on surfaces of response in our work, all j have been set to one in the normalized space. The resulting response surface allows determining the exponent for each cell of the unstructured grid ( Fig. 8 c.). As explained previously, in each experiment the upscaling is done in two or three directions resulting in different K ref ( M, m or t ). Then, two or three values for are obtained per experiment, one for each direction D . In the end, not one but two or three response surfaces of are interpolated and the same amount of omega properties are painted on the unstructured grid. Fig. 8. a) All experiment results, i.e. , are put in the parameter space. b) A simple kriging allows the generation of a surface of response of . c) By projecting cell sizes on the response surface, one is found for every cell of the unstructured grid. Fig. 9. a) One by one, each unstructured cell is filled with points, the number of points depending of the volume of the cell; b) Using the Spectral Turning Bands, permeability is interpolated on these points; c) Using the of the target cell, the permeability of the cell is obtained using power averaging on the points (followed by a back transform from K to K ).

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 . 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 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 resulting from the D-dimensional covariance C . Numerous methods for simulating a onedimensional random field knowing its covariance 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 nonstationary 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. 9 a.). The fine permeability values at these locations are simulated as explained above ( Fig. 9 b.). The global cell value is obtained by averaging the point values using the power law and the estimated for that cell ( Fig. 9 c.). Once again, not one but two or three permeability values are obtained, one for each direction M, m and t .

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.

Permeability simulation
The 3D unstructured grid for this case is composed of 18,250 Voronoï cells with a volume between 25 m 3 and 50.10 3 m 3 ( Fig. 10 a.). 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 25x25x2.5 m and 250x250x25 m for cell sizes varying between 3 and 90 m horizontally and 0.8 to 8 m vertically ( Fig. 10 b.). The porosity has a beta distribution of mean 0.15 and standard deviation 0.05 ( Fig. 10 c.) and the permeability distribution is lognormal with a mean of 100 mD and a standard deviation of 300 mD ( Fig. 10 d.). 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 = = and one for the vertical exponent V ( Fig. 11 ). The two axes of the response surfaces represent the horizontal cell sizes ΔXY (taken as the mean of the cell sizes in X and Y ), and the vertical cell sizes ΔZ . A smooth variation of is observed, from 0 to 1 for H and from -1 to 1 for 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 properties gener- Fig. 10. a) grid for the synthetic case with cell volumes varying with a factor 2.10 3 . 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. Fig. 11. Surfaces of response of horizontal and vertical omega and the corresponding properties on the unstructured grid. Fig. 12. Porosity field on the toy example. Log-permeability will be correlated to this field with a linear coefficient equal to 0.8. ated previously ( Fig. 13 a.). Although both K H and K V properties look similar, the cross-plot on Fig. 13 b. shows that K V is globally inferior to K H . The similarity between these results will be discussed in the conclusion.

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 value given by the response surface). Twenty realizations of porosity and correlated 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 . 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 3 /day and the water (and tracer) injection rate per well is 280 m 3 /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. 16 a.) and the ones of the unstructured grid ( Fig. 16 b.). 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. 16 c.) 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.

Computation times
The use of an unstructured grid, while accounting for the support effect, usually requires to generate the geostatistical simulations on a 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.

Table 1
Computation times for two unstructured grids of different sizes. T 0 is the time to generate the response surfaces and properties on the grids, T 3 is the time to simulate permeability on points using STB and perform power averaging. 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 unstructured 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, 3 = 1 + 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: 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 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.   . 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.

Table 2
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.

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