%~Mouliné par MaN_auto v.0.40.4 (bbb487d5) 2026-04-24 16:48:43
\documentclass[CRMECA,Unicode,biblatex]{cedram}

\addbibresource{CRMECA_Guibert_20260049.bib}

\usepackage{siunitx}
\usepackage{subcaption}
\usepackage{booktabs}

\newcommand{\dif}{\mathop{}\!{\operatorfont{d}}}
\let\ts\textsuperscript

\makeatletter
\renewcommand\p@subfigure{}
\makeatother

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\makeatletter
\def\@setafterauthor{%
  \vglue3mm%
\begingroup\hsize=10.5cm\advance\hsize\abstractmarginL\raggedright
\noindent
   \leftskip\abstractmarginL
  \normalfont\Small
  \@afterauthor\par
\endgroup
\vskip2pt plus 3pt minus 1pt
}
\makeatother

\title{Effects of image resolution and numerical discretization on permeability evaluations}
\alttitle{Effets de la résolution d'image et de la discrétisation numérique sur les évaluations de la perméabilité}

\author{\firstname{Romain} \lastname{Guibert}\CDRorcid{0000-0003-3478-9743}\IsCorresp}
\address{Institut de M\'ecanique des Fluides de Toulouse (IMFT) --- Universit\'e de Toulouse, CNRS-INPT-UPS, Toulouse, France}
\email{romain.guibert@imft.fr}

\author{\firstname{Peter} \lastname{Moonen}\CDRorcid{0000-0002-2551-7669}}
\address{Universit\'e de Pau et Pays de l'Adour, E2S UPPA, CNRS, TotalEnergies, LFCR, Pau, France}
\address{Universit\'e de Pau et Pays de l'Adour, E2S UPPA, CNRS, DMEX, Pau, France}

\author{\firstname{Pierre} \lastname{Horgue}\CDRorcid{0000-0002-2358-8113}}
\address[1]{Institut de M\'ecanique des Fluides de Toulouse (IMFT) --- Universit\'e de Toulouse, CNRS-INPT-UPS, Toulouse, France}

\author{\firstname{Patrice} \lastname{Creux}\CDRorcid{0000-0001-7331-4092}}
\address[2]{Universit\'e de Pau et Pays de l'Adour, E2S UPPA, CNRS, TotalEnergies, LFCR, Pau, France}

\author{\firstname{Franck} \lastname{Plourabou\'e}\CDRorcid{0000-0003-0349-4221}}
\address[1]{Institut de M\'ecanique des Fluides de Toulouse (IMFT) --- Universit\'e de Toulouse, CNRS-INPT-UPS, Toulouse, France}

\author{\firstname{G\'erald} \lastname{Debenest}\CDRorcid{0000-0002-9597-5729}}
\address[1]{Institut de M\'ecanique des Fluides de Toulouse (IMFT) --- Universit\'e de Toulouse, CNRS-INPT-UPS, Toulouse, France}

\thanks{This study was partially funded by Carnot Institute ISIFoR with the ESTIME project. This work was performed using HPC resources from CALMIP (Grant P21016). The Zeiss Xradia 510 used in this work was provided to UAR 3360 DMEX by TotalEnergies}

\keywords{\kwd{Image resolution} \kwd{spatial discretization} \kwd{two-step upscaling} \kwd{permeability}}
\altkeywords{\kwd{Résolution de l’image} \kwd{discrétisation spatiale} \kwd{changement d’échelle à deux étapes} \kwd{perméabilité}}

\begin{abstract}
Digital Rock Physics (DRP) analysis is a widely employed technique for predicting transport parameters from 3D images of core samples. However, the effects of image resolution and spatial discretization of the DRP mesh grid have rarely been systematically studied in detail. To address this issue, we examine a generic sand pack, representing a homogeneous porous medium. This sample was imaged using X-ray micro-tomography at three different spatial resolutions (6, 3, and 1.5~microns/voxel). Permeability is then numerically evaluated by solving the Stokes flow equations using a finite volume method. The processed meshes for converged macroscopic evaluations consist of 105~million to 58~billion cells, necessitating the alternative use of a two-step upscaling method. By employing both methods, this study analyzes the respective influences of image resolution and spatial discretization. Significant effects are observed from both image resolution and spatial discretization, the analysis of which can contribute to identifying optimal strategies for enhancing the accuracy of permeability evaluation.
\end{abstract}

\begin{altabstract}
L'analyse numérique des roches est une technique largement utilisée pour prédire les paramètres de transport à partir d'images 3D d'échantillons. Cependant, les effets de la résolution d'image et de la discrétisation spatiale ont rarement fait l'objet d'études systématiques et approfondies. Pour remédier à cela, nous examinons un paquet de sable générique, représentant un milieu poreux homogène. Cet échantillon a été imagé par micro-tomographie à rayons~X à trois résolutions spatiales différentes (6, 3 et 1,5~micron/voxel). 
La perméabilité est ensuite évaluée numériquement en résolvant les équations de Stokes à l'aide d'une méthode des volumes finis. Les maillages utilisés pour les évaluations macroscopiques comptent entre 105~millions et 58~milliards de cellules, ce qui nécessite le recours à une méthode de changement d’échelle en deux étapes. En utilisant ces deux méthodes, cette étude analyse les influences respectives de la résolution d'image et de la discrétisation spatiale. Des effets significatifs sont observés tant au niveau de la résolution d'image que de la discrétisation spatiale, dont l'analyse peut contribuer à identifier des stratégies optimales pour améliorer la précision de l'évaluation de la perméabilité.
\end{altabstract}

\COI{The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.}

\begin{document}
%\input{CR-pagedemetas}
%\end{document}
\maketitle

\section{Introduction}

Imaging techniques have been widely developed over the past decades, and X-ray Micro-Computed Tomography (MCT) in particular has become very popular~\cite{Wildenschild2013}. MCT allows non-destructive visualization of the micro-structure of many materials in three dimensions~(3D). The technique is used in different fields of application and at various scales. It has become a widely used reference tool in geosciences~\cite{Cnudde2013} where the parallel development of Digital Rock Physics analysis (DRP) \cite{Blunt2013} has provided a large spectrum of numerical methods~\cite{Madonna2012,Andra2013-1,Andra2013-2} to evaluate the transport properties of these materials. MCT provides the detailed morphology of the pore space which provides the basis for flow computations, for example leading to the numerical estimation of permeability.

A major limitation of all imaging methods, including MCT, relies on the intrinsic link between the size of the field of view and the spatial resolution. These two quantities typically differ by three orders of magnitude, since the present size of the camera's image sensor matrix is close to 10\ts{3} pixels in each dimension at the current state of technology and might reach 10\ts{4} in the next decade. This constrained relation between the total field-of-view size and the spatial resolution also has to be compared to the Representative Elementary Volume (REV) length because one would like to evaluate representative quantities from each sample, to obtain statistically meaningful evaluations of transport parameters. When the field of view encompasses the REV length scale, the smallest geometrical features of the available micro-structure or the pore space cannot be less than three orders of magnitude smaller. However, many materials have pore sizes ranging over five or six orders of magnitude, so that many details fall below the image resolution. Conversely, if one attempts to visualize the smallest details, it is then most often impossible for the field of view to reach REV length. In the rare case where the covered range is sufficient, binning techniques can be used to make data more easily handled (see e.g.~\cite{Guan2019} regarding the binning effect on numerical measurements). In general, a compromise must be chosen between resolution and field of view, which is detrimental in many cases. For some materials, an interesting alternative to avoid this aporia is to use so-called statistical methods~\cite{Bazaikin2017}. They consist of simulating the fluid flow through a large number of quasi-representative volumes, imaged at high resolution. This class of methods however cannot guarantee that the threshold of ten voxels per saddle point (the throat region connecting two adjacent pores) is reached and sufficient for the flow to be accurately approximated~\cite{Franc2020}.

Other alternative methods have been recently proposed in the literature~\cite{CHUNG2020106577} using a new numerically efficient approach based on merging voxels rather than using all voxels of the entire three-dimensional CT image. This method locally degrades the initial resolution when clusters of voxels are possible to define. The maximal error found in~\cite{CHUNG2020106577} for permeability evaluation was close to 2{\%} compared to the classical approach without agglomeration with a speedup close to~2.

It is also interesting to mention the recent use of neural networks coupled with imaging techniques to determine porosity, tortuosity, and permeability~\cite{graczyk2020predicting}. Results of the neural network are compared in 2-dimensional images with a conventionally-obtained numerical solution and a good match between the two methods is observed. Recently the same methodology has been applied in 3D~\cite{hong2020rapid} and gives a good match with a speedup close to \num{10000}. In any case, all these new methods are voxel-based methods and, for this reason, cannot avoid the above-mentioned aporia.

A recent attempt to cross-compare voxel-based solver, lattice Boltzmann method, Pore Network Model (PNM), and empirical relations (Kozeny--Carman for instance) on different rock types has been performed~\cite{song2019comparative}. According to them, the results of voxel-based and lattice Boltzmann methods are comparable whereas PNM is less accurate. Correlations overestimate the permeability in all cases. No considerations about image resolution are discussed but standard solvers and pore-scale numerical tools are found to give the most accurate solutions.

A major issue with many popular simulation tools is related to numerical convergence. While the importance of numerical discretization of the pore space is not always emphasized in DRP, it is a well-known central point in computational fluid dynamics. This issue is most often disregarded, especially when a comparative study is performed, on various types of media for example. But this aspect is not to be neglected when one is interested in local phenomena, or precise quantification. Moreover, the effects of discretization on the estimates and observations are a priori directly related to the resolution of the imaged volume.

\pagebreak

The objective of this paper is to quantify the relative impact of numerical discretization on the one hand, and image resolution on the other hand, on permeability estimations of 3D complex porous media. These comparisons are performed on a single sample, chosen to have relatively homogeneous properties. The chosen numerical method is capable of handling very large datasets. The paper is organized as follows: Section~\ref{sec:matmet} describes the 3D porous media sample, its various image acquisitions, their post-processing, and the numerical methods. Section~\ref{sec:rez} describes the results obtained and analyzes the various effects of numerical discretization and image resolution. Conclusions round up the paper.

\section{Materials and methods}\label{sec:matmet}

\subsection{Sample description}

To ensure unbiased results of the current study and avoid potential sample heterogeneity, the investigation is focused on a sand pack. The employed sand originates from Fontainebleau sandstone. It was washed and subsequently dried to retain only the mineral fraction. Nevertheless, no sieving has been performed so a wide range of grain sizes is present. The sand has been retained in a cylindrical polyamide tube of \qty{3.43}{\mm} inner diameter (0.135 inches) with a wall thickness of 76 microns (0.003 inches). Polyamide has the advantage of being almost fully transparent to X-rays.

\subsection{Data acquisition and reconstruction}

We focus on a region of interest (ROI) of \qty{3}{\mm} in diameter and \qty{3}{\mm} in height. A single X-ray tomography of the ROI was performed using a Zeiss Xradia Versa 510 at the DMEX Center for X-ray Imaging (Pau, France). This instrument combines geometrical with optical magnification to yield the desired level of magnification, here targeting a reconstructed voxel size of 1.5 microns. The X-ray source was operated at a \qty{40}{\keV} peak energy. The X-ray beam traverses the sample and the attenuated signal is detected by a 16-bit detector with a 2048$\times$2048 pixel array. A physical filter was employed to limit beam hardening artifacts. In total, 3201 angular projections of the sample were acquired with a \qty{24}{\s} exposure time per projection. These projections were the basis for three reconstructions, targeting 1.5, 3, and 6-micron voxel sizes in the resulting dataset. To that extent binning, i.e.\ the regrouping of pixels, was applied to the raw projections, and the number of projections was reduced accordingly. The resulting datasets are perfectly equivalent to acquisitions of 3201, 1601, or 801 projections of respectively 2048\ts{2}, 1024\ts{2}, or 512\ts{2} pixels per projection. By artificially creating these three datasets from the same collection of raw data, we ensure that we cover each time exactly the same region of interest and that any possible inter-scan variation due to source instability or thermal fluctuations is avoided. The proposed method is quite representative of the magnification changes just idealized to be able to compare the single effect of resolution onto the same dataset.

A horizontal slice taken from the middle of the full 3D dataset is shown in Figure~\ref{fig:samples} for the three considered image resolutions. In these images, dark gray corresponds to the air phase, light gray to the Fontainebleau grains, and nearly white to dense impurities in the sand. Cross-comparison between the different images clearly illustrates that an identical region of interest is visualized at different levels of detail. The 6- and 3-micron resolution images (respectively Figure~\ref{fig:samples}\eqref{fig:samples_a} and Figure~\ref{fig:samples}\eqref{fig:samples_b}) appear blurred compared to the 1.5-micron image (Figure~\ref{fig:samples}\eqref{fig:samples_c}) for which the selected image resolution is appropriate.

\begin{figure}
\begin{subfigure}{.25\textwidth}
\includegraphics[width=\textwidth]{tiff_4_s_250}
\caption{}\label{fig:samples_a}
\end{subfigure}
\hfill
\begin{subfigure}{.25\textwidth}
\includegraphics[width=\textwidth]{tiff_2_s_501}
\caption{}\label{fig:samples_b}
\end{subfigure}
\hfill
\begin{subfigure}{.25\textwidth}
\includegraphics[width=\textwidth]{tiff_1_s_1002}
\caption{}\label{fig:samples_c}
\end{subfigure}
\caption{Middle slice of the sample at different resolutions: \eqref{fig:samples_a}~6~microns, \eqref{fig:samples_b}~3~microns, and \eqref{fig:samples_c}~1.5~microns.}\label{fig:samples}
\end{figure}

\subsection{Image processing}\label{subsec:image-processing}

Permeability calculations are performed on the largest cube, positioned at the center of the visualized volume. Thus the considered volumes in the following are respectively 345\ts{3}, 690\ts{3}, and 1380\ts{3} voxels. The calculations require evaluating the fluid flow through pore space. The latter is extracted from the reconstructed datasets using the open-source software Fiji~\cite{Schindelin2012}. The popular Otsu self-threshold segmentation algorithm~\cite{Otsu1979} has been used to ensure fully reproducible results. Then, isolated, hence inaccessible parts of the pore space are removed using the open-source library CImg~\cite{cimg}, keeping only the largest connected void space. Figure~\ref{fig:binarized} shows the horizontal slice through the middle of the post-processed sample, for the three considered resolutions. We previously discussed in~\cite{en15207796} the anecdotal effect of image filtering for MCT datasets, regardless of whether they are obtained from synchrotron or laboratory sources.

\begin{figure}[b]
\begin{subfigure}{.25\textwidth}
\includegraphics[width=\textwidth]{veryCoarse_345_main_s_173}
\caption{}\label{fig:binarized_a}
\end{subfigure}
\hfill
\begin{subfigure}{.25\textwidth}
\includegraphics[width=\textwidth]{resCoarse_690_main_s_345}
\caption{}\label{fig:binarized_b}
\end{subfigure}
\hfill
\begin{subfigure}{.25\textwidth}
\includegraphics[width=\textwidth]{resFine_1380_main_s_690}
\caption{}\label{fig:binarized_c}
\end{subfigure}
\caption{Middle slice of segmented volume at different resolutions: \eqref{fig:binarized_a}~6~microns, \eqref{fig:binarized_b}~3~microns, and \eqref{fig:binarized_c}~1.5~microns. Grains are shown in black, and pores in white.}\label{fig:binarized}
\end{figure}

Even if pore-space micro-structures are easily recognizable between the three images, one can observe that the porosity of the coarser resolution dataset (Figure~\ref{fig:binarized}\eqref{fig:binarized_a}) is appreciably higher than for the two other datasets. Whilst tailored image processing techniques could certainly improve the segmentation quality of the pore space, especially for low-resolution scans, a differentiated approach would bias the assessment of the image resolution impact on estimated permeability. Since this would undermine the objective of this study, such more elaborated post-processing has been avoided. Nevertheless, as the 6-micron dataset yields an excessively large porosity, it naturally follows that its permeability assessment would be equally excessive and not comparable to the other datasets. For that reason, the 6-micron image is excluded from the numerical comparative study.

\subsection{Models and methods}

\subsubsection{Simple upscaling}

The usual numerical upscaling approach consists of solving the micro-scale Stokes equations (respectively mass and momentum equations)
\begin{gather}
	\nabla\cdot\mathbf{u}=0, \label{eq:stokes-mass}
\\	\mu\Delta\mathbf{u}-\nabla p=\mathbf{0}, \label{eq:stokes-momentum}
\end{gather}
where $\mathbf{u}$ is the micro-scale velocity field, $\mu$ the dynamic viscosity of the fluid, and $p$ is the micro-scale pressure. Then, the permeability $\mathbf{K}$ is deduced, by means of volume averaging, from Darcy's law
\begin{equation}
\mathbf{U}=-\frac{\mathbf{K}}{\mu}\nabla P, \label{eq:darcy}
\end{equation}
where $\mathbf{U}$ is the average fluid velocity (i.e.\ the filtration velocity), $\mathbf{K}$ is the permeability tensor, and $\nabla P$ is the macro-scale pressure gradient. Concerning the boundary conditions, we use the common approach (also previously used in~\cite{Guibert2015,Guibert2016a} consisting in imposing pressure boundary conditions in the direction to be characterized, as well as no-slip boundary conditions along the sidewalls solid surface as in permeameter experiments.
At the Darcy scale, the no-slip condition is represented by a zero-flow boundary, enforced by setting the normal pressure gradient to zero at the walls.
Permeability is numerically evaluated by measuring pressure differences in each direction and the filtration velocity following
\begin{equation}
\mathbf{U} = \bigl\langle \mathbf{u}(\mathbf{x}) \bigr\rangle = \frac{1}{V} \int_{V_p} \mathbf{u}(\mathbf{x}) \dif V,
\end{equation}
where $V$ is the domain volume and $V_p$ is the pore-space volume.

\subsubsection{Two-step upscaling}

This method was initially applied to pore-scale images in~\cite{Horgue2015a}. Synthetically, the Two-Step Upscaling (TSU) consists of (i)~using a usual micro-scale to Darcy-scale upscaling on sub-volumes of the considered sample, (ii)~upscaling the heterogeneous Darcy-scale permeability field obtained at the previous step. This approach has the advantage of considerably reducing both the memory requirement and the computation time while inducing a relatively low error. This induced error can be related to different factors such as rock type, mesh refinement, and sub-volume splitting. In this study, the method is even more interesting, because it allows the processing of large images, whose mesh has been refined several times, which would be almost impossible to process directly.

The computational resources at our disposal allowed us to refine three times the 3-micron image while the 1.5-micron image could only be refined twice. To have the same number of refinement levels for both images, we additionally generated a binned version from the segmented 1.5-micron image, thus having equivalent dimensions to the 3-micron image. Note image processing is non-commutative: binning a segmented image yields another result than segmenting a binned image.

\subsubsection{Numerical methods}

Concerning the meshing step, the initial method~\cite{Horgue2015a} has evolved to allow automated meshing of the various sub-volumes directly from the 3D image, drastically reducing the resources required for the mesh generation phase.

For the numerical resolutions of Stokes' equations, a regular Cartesian mesh is generated respecting the initial image voxels, and retaining only one percolating fluid region (as stated in Section~\ref{subsec:image-processing}). This voxel-based approach, widely used in the literature~\cite{Andra2013-1,Mostaghimi2013}, is reasonable for low-Reynolds flows, i.e.\ in the Stokes regime. Here, we use a dedicated meshing tool that directly generates the fluid domain in a decomposed form, to avoid the memory bottleneck of the native OpenFOAM mesh tools. For the second step of the TSU method, the mesh is simply formed by one cell per subdomain.

The steady-state equations \eqref{eq:stokes-mass}, \eqref{eq:stokes-momentum} and~\eqref{eq:darcy} are solved using the finite-volume toolbox OpenFOAM~\cite{Weller1998}, with second-order spatial discretization schemes and using the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm for Stokes' equations.

\section{Results and discussions}\label{sec:rez}

\subsection{Post-processed data}\label{subsec:post-processed}

Table~\ref{tab:resolution} reports some characteristics of the considered 3D images and numerical meshes. The porosity is defined as the volume of the void space divided by the total sample volume. The porosity for the highest resolution dataset is about 7\% higher than for the 3-micron dataset. For comparison, the characteristics of the 6-micron dataset are also reported. One notices a 50\% change in porosity with respect to the reference image. It highlights that the differences get smaller while increasing resolution. The 1.5-micron image can be considered representative in view of the objective of the current study. The effective porosity, also called connected porosity, is constituted by the pore space which actually contributes to the fluid flow. For the considered sample, the differences between effective and total porosity are negligible and comparable for both resolutions (about 0.06\%).

\begin{table}[!h]
\caption{Key characteristics of the 3D images and the numerical meshes. The characteristics of the 6-micron image are only given for comparison, but the geometry has not been used for the permeability estimate.}\label{tab:resolution}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{@{}ccccc@{}}
\toprule
\textbf{Resolution (\textmu m)} & \textbf{Dimensions~(voxels)} & \textbf{Porosity (--)} & \textbf{\# Cells Reference} & \textbf{\# Cells TSU} \\
\midrule
6 & 345\ts{3} & 0.4741 & not used & not used \\
3 & 690\ts{3} & 0.3225 & 105 876 105 & 105 830 011 \\
1.5 & 1380\ts{3} & 0.3473 & 912 402 405 & 912 294 473\\
\bottomrule
\end{tabular}
\end{table}

\begin{table}[!h]
\caption{Number of cells in global meshes per refinement level.}\label{tab:meshes}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{@{}cccccc@{}}
\toprule
\textbf{Resolution (\textmu m)} & \textbf{L\textsubscript{$-1$}} & \textbf{L\textsubscript{0}} & \textbf{L\textsubscript{1}} & \textbf{L\textsubscript{2}} & \textbf{L\textsubscript{3}} \\
\midrule
3 & -- & 105~M & 846~M & 6.8~B & 54.2~B \\
1.5 & 113~M & 912~M & 7.3~B & 58.4~B & -- \\
\bottomrule
\end{tabular}
\end{table}

The sub-volumes of the TSU method are of identical size for the two images, which are divided into 216 sub-blocks (6$\times$6$\times$6). The decomposition induces an effective fluid volume loss (by creating new unconnected areas) which remains negligible for both image resolutions, of the order of 0.01\% and 0.04\% respectively for the 1.5- and 3-micron datasets. This aspect, specific to the TSU approach, is limited in this study because the sample is homogeneous.

The refined meshes are obtained from the reference mesh ${\mathrm{L}_{0}}$, by a regular 3-directions splitting of the cells. The sizes of the global meshes to be processed are reported in Table~\ref{tab:meshes}. Note that two meshes are formed by more than 50~B cells, which would be very difficult to process directly even on a cluster. This justifies the use of the TSU approach in this case. The refinement level ${\mathrm{L}_{-1}}$ corresponds to the image obtained with binning, which however allows to keep the porosity (with a deviation of about 0.15\%).


Figure~\ref{fig:poro} reports the numerical porosities per slice and per direction. The variations of the porosities are similar at both resolutions. A specific behavior can be observed on the edges linked to the imaging, the shape of the sample, and the packing effects inherent to the constitution of the sample. More specifically, the anisotropy obtained is linked to the average oblong shape of the sand grains which will preferentially position themselves along a horizontal axis during their fall (positioning in the cell) for the equilibrium of grains in the stack. The porosities at the edges are relatively lower than at the center of the sample, which results in a standard deviation of 2--3\% in the transverse directions and more than 5\% in the vertical direction.

\begin{figure}[!h]
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=5cm]{coarsePoroDir}
\caption{}\label{fig:poro_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=5cm]{finePoroDir}
\caption{}\label{fig:poro_b}
\end{subfigure}
\caption{Porosities by direction for \eqref{fig:poro_a}~3-micron resolution volume and \eqref{fig:poro_b}~1.5-micron resolution volume. The vertical direction corresponds to the z-axis.}\label{fig:poro}
\end{figure}

Figure~\ref{fig:diff-images} illustrates the difference between the 3D images obtained with the 1.5- and 3-micron resolutions. These differences appear at the fluid/solid interface, where the gray level varies the most and where thresholding has the largest effect. There is also possibly a centering defect from one image to another because the differences appear more on one side than on the other. However, the objects detected within the packing remain similar, and the differences only occur on their edges, which is the expected effect of the resolution.

\begin{figure}[!h]
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[width=.52\textwidth]{difference_690}
\caption{}\label{fig:diff-images_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[width=.52\textwidth]{difference_690_zoom}
\caption{}\label{fig:diff-images_b}
\end{subfigure}
\caption{Visualization of difference between the 3- and 1.5-micron resolutions: \eqref{fig:diff-images_a}~middle slice of the volume and \eqref{fig:diff-images_b}~zoom in the center region. The volume resolved to 3 microns is interpolated to the volume of the same size as the most resolved one, then the latter is subtracted.}\label{fig:diff-images}
\end{figure}

\subsection{Sub-volume refinement}

In a preliminary step to the full sample processing, and similarly to the mesh convergence performed on a sub-volume in~\cite{Soulaine2016}, we performed simulations on a sub-volume of 100\ts{3} voxels extracted from the most resolved volume. The objective is to illustrate the influence of refinement on the velocity field distribution, and thus the permeability. Figure~\ref{fig:distrib-small-sample} shows the probability density functions obtained from the components of the velocity fields, for a mesh corresponding to the initial resolution (Figure~\ref{fig:distrib-small-sample}\eqref{fig:distrib-small-sample_a}, about \num{400000} cells) a twice refined one (Figure~\ref{fig:distrib-small-sample}\eqref{fig:distrib-small-sample_b}, about 25~M cells). These curves are generated from the velocity fields obtained when resolving the flow in the z direction, which explains the strong di-symmetry of the values represented. Although the general shape of the distributions remains close, two main differences may be identified: (i)~the lowest velocities, which correspond to the fluid/solid interface, are better captured by refining, (ii)~the highest velocities, which correspond to the extreme confinement points (at the level of the saddles points) are also more finely captured. Moreover, a non-negligible shift is observed in the captured maxima.

\begin{figure}[!h]
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=5cm]{100_L0_Unorm}
\caption{}\label{fig:distrib-small-sample_a}
\end{subfigure}
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=5cm]{100_L3_Unorm}
\caption{}\label{fig:distrib-small-sample_b}
\end{subfigure}
\caption{Probability density functions of velocity fields on a sub-volume of 100\ts{3} voxels, for two levels of refinement: \eqref{fig:distrib-small-sample_a}~initial voxel size (${\mathrm{L}_{0}}$) and \eqref{fig:distrib-small-sample_b}~mesh refined two times (${\mathrm{L}_{2}}$).}\label{fig:distrib-small-sample}
\end{figure}

\subsection{Permeability estimation}

Although the method principle is straightforward, simulations on such huge meshes (see Table~\ref{tab:meshes}) remain challenging for technical reasons and thus rather seldom in the literature. We used the Calmip cluster, which is composed of 360~nodes. Each node is itself formed of an Intel\textregistered\ Skylake \qty{2.3}{\GHz} 36-core processor, associated with \qty{192}{\giga\byte} RAM\@. On each sub-volume, depending on the mesh considered, we used between 4~and~20 nodes. The raw permeability data, over the 216~sub-volumes, are reported in Figure~\ref{fig:histo-K}.

\begin{figure}
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=5cm]{coarseDistribL0}
\caption{}\label{fig:histo-K_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=5cm]{fineDistribL0}
\caption{}\label{fig:histo-K_b}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=5cm]{coarseDistribL1}
\caption{}\label{fig:histo-K_c}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=5cm]{fineDistribL1}
\caption{}\label{fig:histo-K_d}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=5cm]{coarseDistribL2}
\caption{}\label{fig:histo-K_e}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=5cm]{fineDistribL2}
\caption{}\label{fig:histo-K_f}
\end{subfigure}
\caption{Histograms and probability density functions of sub-volume permeabilities. Columns correspond to sample resolution: \eqref{fig:histo-K_a},~\eqref{fig:histo-K_c},~\eqref{fig:histo-K_e}~are related to the 3-micron resolution sample and \eqref{fig:histo-K_b},~\eqref{fig:histo-K_d},~\eqref{fig:histo-K_f} to the 1.5-micron resolution sample. Rows correspond to refinement levels: \eqref{fig:histo-K_a}--\eqref{fig:histo-K_b}~initial image resolution (${\mathrm{L}_{0}}$), \eqref{fig:histo-K_c}--\eqref{fig:histo-K_d}~volumes refined once~(${\mathrm{L}_{1}}$) and \eqref{fig:histo-K_e}--\eqref{fig:histo-K_f}~volumes refined twice (${\mathrm{L}_{2}}$).}\label{fig:histo-K}
\end{figure}

The high variability of permeability observed in the raw data highlights that the sub-volumes cannot be considered representative of the full sample. The range of the permeability distribution is directly related to the chosen initial partitioning. In terms of anisotropy, one can observe for every configuration that the vertical permeability is lower than the transverse ones, even if it remains relatively close. Moreover, the refinement decreases the average permeability for both resolutions, as well as reduces the range of sub-volume permeabilities. We can also notice that the decrease of the estimated permeabilities range related to the refinement is less important in the case of the most resolved initial image.

\begin{figure}
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=4.7cm]{hist-coarse}
\caption{}\label{fig:hist-dm_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=4.7cm]{hist-fine}
\caption{}\label{fig:hist-dm_b}
\end{subfigure}
\caption{Histograms of distance maps evaluate from largest component 3D images with \eqref{fig:hist-dm_a}~3-micron resolution and \eqref{fig:hist-dm_b}~1.5-micron resolution.}\label{fig:hist-dm}
\end{figure}

The permeabilities of the full sample for the various configurations, obtained using the second step of the TSU method, are reported in Figures~\ref{fig:K-nat} and~\ref{fig:K-norm}, illustrating the same data but presented in a complementary way. Figure~\ref{fig:K-nat} shows the physical quantities measured versus mesh sizes for each distinct initial image resolution and refinement level. By comparing the results for the two images, it seems that the lowest resolution tends to underestimate the permeability (such as the porosity, see Section~\ref{subsec:post-processed}). The sample anisotropy is observed for various configurations with the permeability in the vertical direction always lower, which confirms a preferential packing effect due to gravity. It is however noticeable that permeability measurements for the largest mesh size (3 microns) are similar for both images, i.e.\ the 3-micron original image and the image constructed by binning from the 1.5-micron image. Despite this similarity, the porosities between the two configurations are different (see Section~\ref{subsec:post-processed}) and one can observe a clear inflection compared to the natural mesh convergence.

Choosing the permeability at the equivalent mesh size to the resolution (without binning) as a reference, we observe a decrease of 21 and 32\% (respectively for the 1.5- and 3-micron mesh size) of the permeability for improved measurements(on two levels). At the same time, the difference in measurements related to the resolution is of the order of 10--12\%.

The slope inflection for the lowest resolutions is explained by a degraded segmentation more hazardous because of a widening of the gray gradients between its pore space and the quartz. Similarly, it can be assumed that the contact between grains is marred, which in turn alters the identification of certain pore thresholds, limiting the connectivity and therefore the permeability.

\begin{figure}
\includegraphics[height=6cm]{k}
\caption{Permeabilities obtained using two-step upscaling function of mesh discretization for the two sample resolutions.}\label{fig:K-nat}
\end{figure}

Figure~\ref{fig:K-norm} shows dimensionless permeability estimation, i.e.\ normalized by the reference permeability value. This reference value $K^0_{\mathrm{ref}}$ is obtained from the direct upscaling method applied on each original mesh to evaluate the difference induced by the TSU method, required for refined meshed. As we previously observed in~\cite{Horgue2015a}, the two-step upscaling induces a permeability reduction of about 3\% for the two images (visible in the figure), which remains relatively small compared to the effect of the other parameters. Results show that the permeability differences between the successive refinement steps are smaller for the finest resolution. This can be explained because a finer resolution tends to reduce the number of pore throats captured by a few cells, typically less than five. As observed in previous work~\cite{Franc2020}, the flow modeling, in this case, is poorly approximated and each refinement step tends to greatly improve the pressure drop estimation. There are no observable links between anisotropy and refinement with similar behaviors for all directions. The trends that emerge from these curves suggest that it would be possible to evaluate the real value once the calculation has converged, by extrapolation.

These numerical results highlight the need to perform some mesh refinement whatever the initial resolution considered to substantially improve the permeability estimation (between 20 to 40\% in our case). This is directly related to the mesh convergence which is generally not performed in this type of simulation and which is however necessary since pore throats may be captured by a few voxels only.

\begin{figure}
\includegraphics[height=6cm]{ktsunorm}
\caption{Semi-log representation of two-step upscaling permeabilities values normalized with the classical approach values at initial image resolution, the function of dimensionless ratio mesh size over voxel size.}\label{fig:K-norm}
\end{figure}

\section{Conclusion}

In this study, a permeability upscaling method has been applied to a sample imaged at several resolutions so as to evaluate and compare the respective influence of mesh size and image spatial resolution on its estimation. The two-step upscaling method has also been used to allow the computation on large meshes, which remains challenging (the largest mesh processed is made of more than 50~billion cells). Although both aspects were expected to have a significant effect on the permeability estimate, a larger effect of the spatial discretization compared to image resolution has been found. Even if the image resolution effect is partly specific to the imaging processing workflow considered in this work, it permits providing some orders of magnitude. These observations allow the set up of numerical simulations depending on the expected objectives. A mere permeability estimator (i.e.\ where one looks for an order of magnitude) can be obtained from a simulation with a Cartesian mesh directly built from the image voxels. On the contrary, if one considers the modeling of transport properties at the micro-scale such as dispersion, this would require several mesh refinements in order to satisfactorily improve the velocity distributions and the pressure drops.

\printCOI

\printbibliography

\end{document}