Geneviève Ndour
email: genevieve.ndour@cea.fr

Tiago Trevisan Jost

Anca Molnos
email: anca.molnos@cea.fr

Yves Durand
email: yves.durand@cea.fr

Arnaud Tisserand
email: arnaud.tisserand@univ-ubs.fr

Evaluation of Approximate Operators Case Study: Sobel Filter Application Executed On An Approximate RISC-V Platform

Keywords: Embedded systems, approximate computing, operators, energy

Approximate computing techniques have been proposed to decrease the energy consumption or improve the performance of computation. Several designs of approximate arithmetic operators are proposed in the literature, suggesting that the operators' energy can be reduced up to 46%. This paper puts existing work in perspective by studying the impact of these operators in a processor targeting the embedded systems. We augment a RISC-V processor with approximate addition and multiplication that support a variable bitwidth in the range of 1 to 16 bits. We investigate the error-energy trade-off for a Sobel filter application. The results indicate that, for an acceptable output image degradation, 7.5% of the energy is saved. We conclude that the benefits of using approximate operators in embedded systems processors, have yet to be proved.

INTRODUCTION

Energy efficiency is a major concern in computing systems. Some emerging applications are error-tolerant, e.g., image processing, object recognition, data mining. Approximate computing is a new field that investigates the trade-offs between application output error and energy or performance [START_REF] Mittal | A survey of techniques for approximate computing[END_REF]. In this context, many approximate operators are proposed. The evaluation of the individual operators indicates large gains in terms of energy consumption.

Embedded systems are frequently constraint in terms of energy. An existing way to exploit the approximate operators work in embedded systems is to integrate the operators themselves or their design principles within hardware accelerators [START_REF] Schlachter | Design of energyefficient discrete cosine transform using pruned arithmetic circuits[END_REF]. In this paper we address a different question, namely the impact of approximate operators in an embedded system processor. For this an applicationlevel evaluation is needed.

Existing application level evaluations of approximate operators are limited to the floating point number representation [START_REF] Chandrasekharan | ProACt: A Processor for High Performance On-demand Approximate Computing[END_REF][START_REF] Sampson | EnerJ: Approximate data types for safe and general low-power computation[END_REF]. However embedded systems often makes use of integer operators. Subsequently we augment a RISC-V processor [7] with approximate addition and multiplication and we investigate the errors due to Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from permissions@acm.org. SAMOS XVIII, July 2018, Samos Island-Grece © 2016 Association for Computing Machinery. ACM ISBN 123-4567-24-567/08/06. . . $15.00 https://doi.org/10.475/123_4 approximations on a Sobel filter case study. The instruction breakdown of the Sobel filter indicates that the application contains 93.2% of non-approximable operations, such us memory accesses, and branches. We perform several optimizations to increase the number of approximable operations. The error-energy trade-off for the Sobel filter application, is performed on a bit range relevant for the image processing applications, i.e., 1 to 16 bits. The error metrics of interest for the Sobel filter are: the root mean square error (RMSE) [START_REF] Sivakumar | Comparative study of speckle noise reduction of ultrasound b-scan images in matrix laboratory environment[END_REF], the peak signal to noise ratio (PSNR) [START_REF] Sivakumar | Comparative study of speckle noise reduction of ultrasound b-scan images in matrix laboratory environment[END_REF] and the structural similarity (SSIM) [START_REF] Wang | Structural similarity based image quality assessment[END_REF]. Finally, we observe that, for an acceptable output image degradation, 0.95 SSIM, only 7.5% of the energy is saved.

The rest paper is organized as follows. Section 2 presents the related work, Section 3 describes the target architecture and the energy model, Section 4 presents the experimental results on the Sobel filter program and Section 5 concludes the paper.

RELATED WORK

Several approaches have been proposed in approximate computing at different layers to reduce the computation costs (e.g., energy). For example, instruction set architecture (ISA) extensions with approximate operations aim to reduce energy consumption [START_REF] Hadi Esmaeilzadeh | Architecture support for disciplined approximate programming[END_REF]. In [START_REF] Chandrasekharan | ProACt: A Processor for High Performance On-demand Approximate Computing[END_REF][START_REF] Sampson | EnerJ: Approximate data types for safe and general low-power computation[END_REF] the authors propose to utilize 32-bit approximate floating point operators that allow the scaling of the operands bit-width. These studies present application-level error-energy trageoffs, however not for integer operators.

Several integer approximate operators have been also designed. In [START_REF] Jahier Pagliari | A methodology for the design of dynamic accuracy operators by runtime back bias[END_REF], Pagliari et al. propose a 16-bit multiplier using supply voltage scaling to reduce the bit-width of the operator. As a result energy consumption is reduced by up to 39%. In [START_REF] Andrew | Accuracy-configurable adder for approximate arithmetic designs[END_REF], Kahng et al. propose an accuracy-configurable approximate (ACA) adder which reduces energy consumption up to 30%. In [START_REF] Khaing Yin Kyaw | Low-power highspeed multiplier for error-tolerant application[END_REF], Kyaw et al. propose an error tolerant multiplier (ETM) that composed of two parts: a most significant bit (MSB) part for high accuracy and a least significant bits (LSB) part that allows certain level of error. The 12-bit ETM reduces energy consumption by up to 46%. We can observe that the energy savings in the above individual operators can achieve more than 30%. However, these operators have not been evaluated on a processor, at the application-level, to measure the global energy gains. The closed related work is the study of approximate floating point on the Sobel filter [START_REF] Chandrasekharan | ProACt: A Processor for High Performance On-demand Approximate Computing[END_REF]. The results indicate energy gains with little impact in output quality. Nevertheless, this study is limited to precisions above 20 bits, and hence the most interesting bitwidths for image processing, i.e., under 16 bits, are not explored. To the best of our knowledge, we are the first to integrate approximate integer operator in a RISC-V processor and to evaluate their impact on an image processing application, namely a Sobel filter.

TARGET ARCHITECTURE

The target architecture is a reduced instruction set computer (RISC) core, which is often used in low-power embedded systems. We employ the RISC-V 32-bit architecture template and supplement it with approximate operations. In what follows, the extension of RISC-V architecture and the energy consumption model are presented.

RISC-V architecture extended with approximate operations

For ease of explanation, we consider a classical 5-stage RISC-V pipeline. The pipeline stages are: IF (instruction fetch) to bring the instruction from the memory; ID (instruction decode) to decode the instruction; EXE (execution) to execute the instruction; MEM (memory) for data transfer instructions; WB (write back) to store the output data into the destination register. Note that our investigation is valid for other possible partitions of a RISC processor into pipeline stages.

We extend this RISC-V architecture by adding in the pipeline a unit to execute the approximate operators, denoted as a.EXE, as presented in Figure 1. a.EXE implements addition and multiplication with a configurable degree of approximation. In this work, approximation is implemented by executing the operation on a reduced bit-width, i.e., only on a configured number of most significant bits [START_REF] Jahier Pagliari | A methodology for the design of dynamic accuracy operators by runtime back bias[END_REF]. A separate instruction was added to configure the bit-width for a.EXE. All other RISC-V instructions are denoted as accurate instructions.

Energy model

To estimate the energy consumption of an application, we consider an energy model per instruction. This model is constructed starting with measurements obtained from a test-chip that includes a processor core with a similar pipeline as the RISC-V one and 256KB of memory, split into 64KB of instruction memory and 192KB of data memory. The energy is an average over multiple executions with random input data, and we assume that the energy of an instruction does not depend on the preceding instructions. In the rest of this paper, we present the energy values relative to multiplication (1) To simplify the model, we split the RISC-V instructions into classes with similar energy consumption. For the accurate instructions the classes are: additions (add) that include additions and subtractions; multiplications (mul); logical instructions (logic), e.g., xor, and; branch instructions (br), e.g., jump, branch on greater than (bgt);

where e C PU bas e is the energy consumed on the stages IF, ID, MEM and WB in Figure 1, e C PU a .EX E (bw ) is the energy consumed to execute an approximate instruction. e C PU a .EX E (bw ) is obtained by fitting energy values from literature, in our case [START_REF] Jahier Pagliari | A methodology for the design of dynamic accuracy operators by runtime back bias[END_REF], such that e C PU a .EX E (bw ) is equal to 0 for 0 bit and for the maximum bitwidth, e C PU a .EX E (bw ) is calculated with the following formula:

e C PU a .EX E (maxbw ) = e C PU EX E × (1 + f aovh ) (3)
where f aovh > 0 is the overhead energy of implementing the approximate operators due to the partitioning of the operators into various domains with different supply vector threshold, that introduces area overheads as we can see in [START_REF] Jahier Pagliari | A methodology for the design of dynamic accuracy operators by runtime back bias[END_REF]. In our case, f aovh = 0.1. Figure 2 presents the relative energy values of the approximate instructions (a.add and a.mul), depending on the bit-width.

EXPERIMENTAL CASE STUDY: THE SOBEL FILTER

The Sobel filter is a kernel used in image processing for edge detection. Edge detection is an image processing technique that identifies object boundaries in an image. The input is a RGB image and the output is a gray-scale image in which the edges are emphasized. The image gradient of each pixel is calculated by convolving the image with a pair of filters (horizontal and vertical filters), that are 3x3 matrix.

To perform experiments, the RISC-V spike simulator [START_REF]Spike simulator[END_REF] is augmented with the capabilities to count instructions and energy per class. The aim is to determine the fraction of the instructions that can be aproximated and investigate the error-energy tradeoff, as presented in the following.

Instructions breakdown

We utilize the Sobel filter source code provided by the Axbench benchmark [START_REF] Yazdanbakhsh | AxBench: A Benchmark Suite for Approximate Computing Across the System Stack[END_REF]. However this implementation of the Sobel algorithm is in floating point and it is not optimised for the execution on a embedded processor. Consequently, we converted the source code from floating point to integer and we applied several optimizations typical for the low power domain. The optimisations are stacked-up one after the other, as follows:

• V1.: 2D vector for pixel storage is replaced by 1D vector to reduce the performance costs due to address calculations. • V2.: V1 + the redundant code lines are removed to avoid to compute several times the same value, because the 3 components of a pixel have the same value in gray-scale image. • V3.: V2 + the loops for convolution computations are unrolled to reduce the instructions introduced by the the loop indexing. • V4.: V3 + remove filtering with null values. To apply the filters on the pixels that are at the border of the image, in the original version, the image is padded with zeros. In this optimisation we rewrite the filtering loops to remove this padding. In Figure 3 we present the effect of these optimisations in terms of instructions breakdown. The percentage of approximate instructions is equal to 6.8% and 17% in the original version and the most optimized one, V4, respectively. Despite optimizations, most instructions are still not approximable (83% in the optimized one). These instructions consist of accurate arithmetic for address calculation, data memory access and branch instructions. The number of load instructions are generally large in comparison to the others. The high percentage of the load instructions is due to the fact that the pixels of the image are stored into the data memory and to compute the gradient of one pixel, the 3x3 matrix elements of the two filters and ones of the region of the pixel are loaded from the data memory.

In the next section, the most optimized version, V4, is utilized in the experiments.

Error-energy trade-off

To evaluate the errors induced to the Sobel filter by the reduce bit-width, we utilize three of the most common metrics proposed in image processing domain: the root mean square error (RMSE) [START_REF] Sivakumar | Comparative study of speckle noise reduction of ultrasound b-scan images in matrix laboratory environment[END_REF], the peak signal to noise ratio (PSNR) [START_REF] Sivakumar | Comparative study of speckle noise reduction of ultrasound b-scan images in matrix laboratory environment[END_REF] and the structural similarity (SSIM) [START_REF] Wang | Structural similarity based image quality assessment[END_REF]. The RMSE and the PSNR estimate the absolute errors between pixels of two compared images but are less correlated to human perception of image quality compared to the SSIM. To calculate the errors we consider as reference image the output of the original floating point Sobel filter. The experiments are performed on 100 input images.

Figure 4 presents the mean and the standard deviation of the RMSE and PSNR metrics with various bit-widths; the lower the RMSE, the higher the output quality is; the higher the PSNR, the Figure 5, evaluates the similarity between the reference and the approximate images. The results indicate that at least 8 bits are required to have an output image similar to the reference image.

Figure 6 shows that the energy gains vary according to the quality of the output and to have an acceptable output, the energy saved is equal to 7.5%. The results indicate that, even after code optimization, the energy gains are not as large as the ones obtained evaluating the individual operators separately, as we can see in Table 2. The reduction of energy savings is due to the fact that the application includes more than 80% of non-approximate instructions.

CONCLUSION AND DISCUSSION

This paper presents an evaluation of common approximate operations, i.e., addition and multiplication, in a RISC-V processor, on a Sobel filter edge detection use-case. The operations are of integer type and the approximation therein is implemented by varying the bit-width of the operations. The Sobel output quality is evaluated with several metrics and the energy consumption is estimated. We employ an energy per instruction model derived from measurements on a test-chip.

The results indicate that even after optimizations, the fraction of the approximable instructions in the Sobel filter remains low, i.e., 17%. The memory accesses, the arithmetic instructions for computing addresses, and the branches, account for the largest fraction of the Sobel filter instructions and are not approximated. Furthermore, the Sobel filter, does not require more than 16 bits to execute with no errors when compared to a floating point implementation. This reduces the possible range of bit-width scaling, and hence of energy scaling. Note, this is the case for more image processing filters. Finally, the energy consumption of approximate addition and multiplications contains a part that does not scale down when the bit-width is reduced. These explain why the maximum energy savings are less than 10% for the Sobel filter.

To exploit approximate operations in embedded processors several directions for investigation are in sight. One is to extend the bit-width reduction also when accessing the data memory, i.e., load/store only the MSB bits. Note the other dominant instructions, i.e., branches, address computation, cannot be easily approximated without severe errors. Second, the study should be extended beyond 16 bits for embedded applications that require such a bitwidth. Third, larger approximate operators, e.g., fused, have a higher chance of bringing energy gains without a significant impact in the quality of result.

1

 1

Figure 1 :

 1 Figure 1: RISC-V Architecture

  . The energy (e t ot al ) consumed by an instruction is the sum of the energy consumed by the instruction memory (e MemI nst ), the energy consumed by the CPU (e C PU ) and the energy consumed by the data memory (e MemDat a ), as follows: e t ot al = e MemI nst + e C PU + e MemDat a

Figure 2 :

 2 Figure 2: Relative energy values of a.add and a.mul

Figure 3 :Figure 4 :

 34 Figure 3: Sobel filter instructions

Figure 5 :Figure 6 :

 56 Figure 5: SSIM between the reference image and the approximate image

Table 1 :

 1 Relative energy values of the accurate instructions The test chip allows separate measurements for the energy of the memory blocks and CPU. The total energy, e tot al , is presented per class on Table1. The approximate additions (a.add) and multiplications (a.mul) are considered to be each in a class of their own, to gain insights into their utilization. For these instructions, the e C PU depends on the bit-width (bw), as follows:

	Classes of instructions relative energy values
	add	0.69
	mul	1
	logic	0.67
	br	1.56
	st	1.27
	ld	0.9
	store (st) and load (ld).

e C PU (bw ) = e C PU bas e + e C PU a .EX E (bw ),

Table 2 :

 2 Energy savings

	Approximate operators/programs Energy savings
	Adder 16-bit [3]	30%
	Multiplier 12-bit [4]	46%
	Multiplier 16-bit [6]	39%
	Sobel filter	less than 10%

ACKNOWLEDGEMENTS

This work was partially supported by the French Agence Nationale de la Recherche, under Grant Agreement ANR-15-CE25-0015, project ARTEFACT. Furthermore, we would like to thank our colleagues from the Design and Architecture Lab for providing us with the energy values.


Philippe Dagaut

Yuri Bedjanian

Guillaume Dayma

Fabrice Foucher

Benoît Grosselin

Manolis Romanias

Roya Shahla

EMISSION OF CARBONYL AND POLYAROMATIC HYDROCARBON POLLUTANTS FROM THE COMBUSTION OF LIQUID FUELS: IMPACT OF BIOFUEL BLENDING

The combustion of conventional fuels (Diesel and Jet A-1) with 10-20% vol. oxygenated biofuels (ethanol, 1-butanol, methyl octanoate, rapeseed oil methyl ester, diethyl carbonate, tri(propylene glycol)methyl ether, i.e., CH 3 (OC 3 H 6 ) 3 OH, and 2,5-dimethylfuran) and a synthetic paraffinic kerosene was studied. The experiments were performed using an atmospheric pressure laboratory premixed flame and a four-cylinder four-stroke Diesel engine operating at 1500 rpm. Soot samples from kerosene blends were collected above a premixed flame for analysis. Polyaromatic hydrocarbons (PAHs) were extracted from the soot samples. After fractioning, they were analyzed by highpressure liquid chromatography (HPLC) with UV and fluorescence detectors. C 1 to C 8 carbonyl compounds were collected at the Diesel engine exhaust on 2,4dinitrophenylhydrazine coated cartridges (DNPH) and analyzed by HPLC with UV detection. The data indicated that blending conventional fuels with biofuels has a significant impact on the emission of both carbonyl compounds and PAHs adsorbed on soot. The global concentration of 18 PAHs (1-methylnaphthalene, 2-methyl-naphthalene, and the 16 US priority EPA PAHs) on soot was considerably lowered using oxygenated fuels, except 2,5-dimethylfuran. Conversely, the total carbonyl emission increased by oxygenated biofuels blending. Among them, ethanol and 1-butanol were found to increase considerably the emissions of carbonyl compounds.

INTRODUCTION

Combustion of conventional fuels such as gasoline, Diesel or jet fuels currently used for transportation significantly contributes to greenhouse emissions and global warming [START_REF] Sawyer | Science based policy for addressing energy and environmental problems[END_REF]. Recently the interest for synthetic and bio-derived fuels [START_REF] Escobar | Biofuels: Environment, technology and food security[END_REF][START_REF]Biofuels flight path[END_REF][START_REF] Mzé Ahmed | The Oxidation of a Coal-to-Liquid Synthetic Jet Fuel: Experimental and Chemical Kinetic Modeling Study[END_REF][START_REF] Dagaut | Combustion of a Gas-to-Liquid-Based Alternative Jet Fuel: Experimental and Detailed Kinetic Modeling[END_REF][START_REF] Dagaut | Experimental and detailed kinetic model for the oxidation of a Gas to Liquid (GtL) jet fuel[END_REF][START_REF] Dagaut | The Combustion of Synthetic Jet Fuels (Gas to Liquid and Coal to Liquid) and Multi-Component Surrogates: Experimental and Modeling Study GT2015-42004[END_REF][START_REF] Mze-Ahmed | Kinetics of Oxidation of a 100% Gas-to-Liquid Synthetic Jet Fuel and a Mixture GtL/1-Hexanol in a Jet-Stirred Reactor: Experimental and Modeling Study[END_REF][START_REF] Dagaut | Experimental and Kinetic Modeling of the Oxidation of Synthetic Jet Fuels and Surrogates[END_REF][START_REF] Dagaut | Combustion of synthetic jet fuels: Naphthenic cut and blend with a gas-to-liquid (GtL) jet fuel[END_REF], which are considered helpful for reducing dependence of both air and ground transportation on petroleum [START_REF] Hermann | Chemical analysis of combustion products from a high-pressure gas turbine combustor rig fueled by jet A1 fuel and a Fischer-Tropsch-based fuel[END_REF][START_REF] Corporan | Emissions characteristics of a turbine engine and research combustor burning a Fischer-Tropsch jet fuel[END_REF][START_REF] Huber | Surrogate mixture model for the thermophysical properties of synthetic aviation fuel S-8: Explicit application of the advanced distillation curve[END_REF], and environment friendly, increased. The Fischer-Tropsch process allows the production of kerosene or Diesel types of fuel from synthesis gas. Nevertheless, exhaust emissions released into the atmosphere lead to major environmental and health concerns. Some polluting constituents such as carbon monoxide, nitrogen oxides, unburned hydrocarbons, and particulate matter [START_REF] Buchholz | Using Carbon-14 Isotope Tracing to Investigate Molecular Structure Effects of the Oxygenate Dibutyl Maleate on Soot Emissions from a DI Diesel Engine[END_REF][START_REF] Ratcliff | Knock Resistance and Fine Particle Emissions for Several Biomass-Derived Oxygenates in a Direct-Injection Spark-Ignition Engine[END_REF] are strictly regulated by governments. However, other exhaust constituents, generally referred to as unregulated pollutants, have significant impact on atmospheric pollution and human health. Carbonyl compounds (CBCs) and polyaromatic hydrocarbons (PAHs) are among the important unregulated pollutants emitted from internal combustion engines and gas turbines. Carbonyl compounds, i.e., aldehydes and ketones, can be produced from incomplete combustion of fuels [START_REF] Ratcliff | Impact of Higher Alcohols Blended in Gasoline on Light-Duty Vehicle Exhaust Emissions[END_REF] and are of significant importance in atmospheric chemistry [START_REF] Seinfeld | Atmospheric Chemistry and Physics: From Air Pollution to Climate Change[END_REF] through their oxidation and interaction with nitrogen oxides which yield other pollutants, e.g., ozone and peroxyacyl nitrates, lowering urban air quality. In addition to their irritative and ecotoxic properties, some carbonyl compounds such as formaldehyde, acetaldehyde and acrolein are identified as probable or known carcinogen and/or mutagen. Polycyclic aromatics originating principally from incomplete combustion of fossil fuels and other organic matter are present on particulates and in the gas phase. PAHs and their derivatives adsorbed on soot have been largely studied due to their potential carcinogenic and/or mutagenic effects, and reactivity in the troposphere. While there is an increasing interest for developing and using biofuels with the aim of reducing dependency to oil, improving engine efficiency and reducing pollutants formation are needed. Oxygenated and synthetic fuels used pure or as additive are potential alternatives to petroleum-based fossil fuels with expected, but not well-assessed, beneficial effects on emissions.

The purpose of this work is to study the impact on pollutants emissions of blending conventional liquid fuels (Jet A-1 and Diesel) with synthetic jet-fuel and oxygenates that can derive from biomass. The pollutants considered are polyaromatic hydrocarbons, ketones and aldehydes. Carbonyls emissions were measured in the exhaust of a Diesel engine running with Diesel and blends. The identification and quantification of PAHs adsorbed on soot were performed using a premixed flame and Jet A-1, easier to handle with a laboratory burner, as reference fuel. Due to the similar composition of Diesel and Jet A-1 fuels, it is expected that there will be no major impact on the interpretation of the present data.

EXPERIMENTAL

For engine experiments, a conventional Diesel fuel (DF) from Caldeo was used as reference fuel. The other fuels were prepared by blending the DF with 10% vol. of the following oxygenated additives: ethanol (EtNOL), 1-butanol (1-BNOL), methyloctanoate (MOC), diethylcarbonate (DEC), tri(propylene glycol) methyl ether (TPGME), rapeseed oil methyl ester (RME), and a synthetic paraffinic kerosene (SPK) from Sasol having a derived cetane number of 23 (ASTM D7668). Derived cetane numbers of Diesel fuel and fuel blends tested in engine were measured using a Herzog Cetane ID 510 Analyzer (ASTM D7668); they are given in Table 1 together with the average Filter Smoke Number (FSN) which was continuously measured using an AVL-450s smoke meter in the experiments. The engine used in this study is a direct-injection 4cylinder 4-stroke Diesel engine with a Continental common-rail injection system. The exhaust gas recirculation rate was set to 32.5%. Experiments were performed at a steady engine operating conditions: BMEP = 0.75 MPa (40% engine load) under a constant speed of 1500 rpm. The main engine characteristics are given in Table 2. Before refilling the engine with a new fuel, the connection tubing and fuel tank were drained out and the fuel filter was replaced. Prior to each experiment, the engine was allowed to warm up for 30 min to ensure the remaining fuel in the system was consumed, and to stabilize the operating parameters. The experiments were repeated at least once. No aftertreatment system was used. Sampling and analysis of CBCs were based on the US EPA TO-11A method [START_REF]Compendium method TO-11A: Determination of formaldehyde in ambient air using adsorbent cartridge followed by high performance liquid chromatography (HPLC) [Active sampling methodology[END_REF]. Ten gas samples were collected at 3 m away from the outlet of the combustion chamber by using an automated sampling device (Figure 1). The exhaust gases were collected on 2,4-dinitrophenylhydrazine (DNPH)-coated silica gel cartridges at a constant sampling flow of 0.5 L/min for 5 min. DNPH reacts with carbonyls to form the corresponding stable dinitrophenylhydrazone derivatives that are quantified. At the aforementioned sampling distance from the engine, the gas temperature is low enough (ca. 60°C) to avoid damaging DNPH or the carbonyls derivatives. In order to minimize errors due to the variation of sampling flow rate, ten samples were collected for each fuel, and the average values were used.

DNPH cartridge

After sampling, the cartridges were slowly eluted with 5.0 mL acetonitrile and the extracts were labeled and stored at 4 °C. The samples containing the carbonyl-DNPH derivatives were then analyzed by high-performance liquid chromatography. HPLC grade acetonitrile (ACN), water and tetrahydrofuran (THF) were employed as mobile phase. The following gradient elution was adopted for ACN/water/THF: 0-12 min at 40:55:5 v/v/v and a flow rate of 1.0 mL/min followed by a linear gradient to 55:40:5 v/v/v until 35 min and a flow rate of 1.5 mL/min, and held at 55:40:5 v/v/v until 50 min with a flow rate of 1.5 mL/min. The sample injection volume and column temperature were 20 µL and 40 °C, respectively. The carbonyls detection wavelength was set to 360 nm. The identification and quantification of CBCs were performed using a DNPH-carbonyl standard containing 13 carbonyl compounds: formaldehyde, acetaldehyde, acrolein, acetone, propionaldehyde, crotonaldehyde, methacrolein, 2-butanone, n-butyraldehyde, benzaldehyde, valeraldehyde, p-tolualdehyde, and hexanal. For premixed flat-flame burner experiments, a conventional aviation jet fuel (Jet A-1) with a derived cetane number of 44.7 (ASTM D7668) was used as reference fuel. Four fuel blends were prepared by adding 20% by volume of the following oxygenated additives to Jet A-1: 1-butanol, methyl octanoate, diethyl carbonate, 2,5-dimethylfuran, and a synthetic paraffinic kerosene from Sasol. A laboratory flat-flame burner was used for the production and deposition of soot samples from premixed flames of liquid fuels (Jet A-1 and fuel blends) in a reproducible way and under well-characterized and controlled combustion conditions [START_REF] Lelièvre | Heterogeneous reaction of ozone with hydrocarbon flame soot[END_REF]. The soot samples were produced under fuel-rich condition (φ ≈ 2.3) and deposited on the outer surface of a Pyrex tube. The collection tube was introduced perpendicular to the flame axis at the height of 4 cm above the head of the burner and was rotated and moved through the flame by operator. The flame was isolated from environment atmosphere with a cylindrical Pyrex tube through which a nitrogen flow was introduced from the small orifices creating a pressure lightly higher than atmospheric pressure. This avoids the air diffusion and helps to have a homogeneous and stabilized flame. The sampling Pyrex tube was thermostated at 45 °C by an internal water circulation in order to prevent surface temperature gradient and ensure a homogeneous distribution of PAHs along the tube (Figure 2).

RESULTS AND DISCUSSION

Engine results

Thirteen carbonyl compounds were identified and quantified in exhaust emissions of a four-cylinder four-strike Diesel engine running at a stabilized operating point with conventional DF and 7 blends (10% by volume of biofuels or SKP to limit the impact on engine operating conditions such as combustion phasing).

The measured total CBCs emissions are shown in Figure 3. The total concentration of carbonyl compounds obtained for the Diesel fuel was close to 1 mg m -3 . Based on our results, blending DF with RME and SPK has no major effect on the emission of carbonyl compounds. However, blending DF with ethanol, 1-butanol, diethylcarbonate, methyloctanoate, and TPGME increases significantly the concentrations carbonyl compounds in exhaust gases. Here, the highest CBCs content was measured for the DF/Ethanol blend (3.87 mg m -3 ). Carbonyl compounds contribution in exhaust emissions is shown in Figure 4. Under the operating conditions applied in this study, acetaldehyde was the most abundant carbonyl compound measured in the exhausts, for all fuels. Acetaldehyde emissions were in the range 0.32-1.72 mg m -3 , which corresponds to 33-45% of the mass of carbonyl compounds in exhaust gases. In addition, the concentration of acetaldehyde in exhaust emissions is higher than formaldehyde for all fuels. Formaldehyde, with a highest contribution of 0.3 mg m -3 , is the second most abundant carbonyl compound found after acetaldehyde. It represents 6-14% of the total mass of carbonyl compounds in exhaust gases. The emissions of individual carbonyl compounds are presented in Figure 5. As shown in Figure 5-(a), acetaldehyde is the major component among carbonyls followed by formaldehyde and acrolein. However, one should note that the measured concentrations of carbonyl compounds are very low (less than 1 ppm). A global uncertainty of ± 25% was estimated. Acetone, propionaldehyde, and crotonaldehyde were also emitted in significant quantities, but their concentrations were below those of the three major compounds cited above. Similar results have been reported by Guarieiro et al. [START_REF] Guarieiro | Emission profile of 18 carbonyl compounds, CO, CO2, and NOx emitted by a diesel engine fuelled with diesel and ternary blends containing diesel, ethanol and biodiesel or vegetable oils[END_REF] for carbonyl compounds emitted by a Diesel engine operating under similar conditions (1800 or 2000 ppm).

Formaldehyde and acetaldehyde emissions increased for all blended fuels. The highest emissions of formaldehyde (0.223 ppm) and acetaldehyde (0.874 ppm) were measured for the DF/DEC blend and the DF/Ethanol blend, respectively.

Regarding the DF/Ethanol blend, formaldehyde, acetaldehyde and acetone were the most abundant carbonyls. Similar results regarding the increase of formaldehyde and acetaldehyde emissions by adding ethanol to Diesel fuel have been reported by Song et al. [START_REF] Song | Carbonyl compound emissions from a heavy-duty diesel engine fueled with diesel fuel and ethanol-diesel blend[END_REF], and Cheung et al. [START_REF] Cheung | Experimental investigation of regulated and unregulated emissions from a diesel engine fueled with ultralow-sulfur diesel fuel blended with ethanol and dodecanol[END_REF]. The partial combustion of the fuel and the presence of ethanol as the main precursor of acetaldehyde (CH 3 CH 2 OH + • OH → CH 3

• CHOH → CH 3 CHO) can explain the high concentrations of formaldehyde and acetaldehyde [START_REF] Agarwal | Biofuels (alcohols and biodiesel) applications as fuels for internal combustion engines[END_REF]. As proposed by Tao et al. [START_REF] Tao | Investigation of the chemical structures of laminar premixed flames fueled by acetaldehyde[END_REF], the degradation of acetaldehyde can result in the formation of CH 3

• CO that can lead to the formation of acetone via:

CH 3 • CO + • CH 3 → CH 3 COCH 3
Concerning the DF/DEC blend, acetaldehyde was again the most abundant carbonyl emitted, followed by formaldehyde. Such important quantity for acetaldehyde is consistent with experimental results obtained from the oxidation of diethyl carbonate in a jet-stirred reactor [START_REF] Shahla | Burning velocities and jet-stirred reactor oxidation of diethyl carbonate[END_REF] where acetaldehyde formation can result from:

(CH 3 CH 2 O) 2 C=O → CH 3 CH 2 O • → CH 3 CHO (CH 3 CH 2 O) 2 C=O → CO 2 + C 2 H 4 + CH 3 CH 2 OH CH 3 CH 2 OH + • OH → CH 3 • CHOH → CH 3 CHO + H • Figure 5-(b)
focuses on other carbonyls (methacrolein, 2-butanone, n-butyraldehyde, benzaldehyde, valeraldehyde, p-tolualdehyde and hexanal) which were present in the exhaust emissions but in much lower concentrations. High emissions of nbutyraldehyde and 2-butanone were observed with the DF/1-Butanol blend. The presence of the decomposition products of 1-butanol in the reaction zone could result in the formation of these two carbonyl compounds:

CH 3 (CH 2 ) 3 OH + X → CH 3 (CH 2 ) 2 • CHOH + XH CH 3 (CH 2 ) 3 OH +X → CH 3 CH • CH 2 CH 2 OH + XH
which can in turn yield n-butyraldehyde and 2butanone. n-Butyraldehyde is formed via:

CH 3 (CH 2 ) 2 • CHOH → CH 3 (CH 2 ) 2 CHO + H
The formation of 2-butanone can derive from the lowtemperature oxidation of 1-butanol triggered by that of the Diesel fuel (with multiple oxygen addition to fuelderived radicals, internal hydrogen transfer, and decomposition, followed by recombination of • CH 3 and CH 3 (CO) • CH 2 ):

CH 3 • CHCH 2 CH 2 OH (+ 2 O 2 + CH 3 ) → CH 3 C(O)CH 2 CH 3 (+ • OH + HCOOH)
The seven carbonyls shown in Figure 5-(b) contribute to 10-28% of the total emissions.

Premixed flame results

The soot particulates were produced and collected using an atmospheric pressure premixed flame of Jet A-1, SPK, and four biofuel-Jet A-1 blends under fuelrich conditions (φ ≈ 2.3). The soluble PAHs fractions were extracted and analyzed by HPLC system using a recently validated analytical methodology [START_REF] Andrade-Eiroa | An alternative to trial and error methodology in solid phase extraction: an original automated solid phase extraction procedure for analysing PAHs and PAHderivatives in soot[END_REF]. Based on our first observations on premixed flame, the sooting tendency of fuels can be assumed proportional to soot collection rates defined as the total mass of soot particulate deposited in time. One should note that this trend cannot be simply estimated from the O/C ratio in the fuel (Table 3) because it depends strongly on the fuel chemical structure [START_REF] Mcenally | Fuel decomposition and hydrocarbon growth processes for oxygenated hydrocarbons: butyl alcohols[END_REF][START_REF] Mcenally | Sooting Tendencies of Oxygenated Hydrocarbons in Laboratory-Scale Flames[END_REF][START_REF] Pepiot-Desjardins | Structural group analysis for soot reduction tendency of oxygenated fuels[END_REF]. Indeed, blends containing 1-butanol or methyl octanoate produce significantly different amounts of soot whereas they have the same O/C ratio, which is in line with earlier findings showing alcohols are more efficient to reduce soot than esters [START_REF] Buchholz | Using Carbon-14 Isotope Tracing to Investigate Molecular Structure Effects of the Oxygenate Dibutyl Maleate on Soot Emissions from a DI Diesel Engine[END_REF]. The relative sooting tendencies derived here from this assumption were: 1 for Jet A-1, 0.33 for SPK, 0.37 for Jet A-1/1-butanol, 0.50 for Jet A-1/DEC, 1.17 for Jet A-1/Methyloctanoate, and 1.77 for Jet A-1/2,5dimethylfuran (2,5-DMF). Table 4 shows the measured absolute concentrations of the 16 US EPA priority PAHs present on soot sample. The total PAHs concentration varied from 250 to 4700 ng mg -1 . Conventional Jet A-1 and SPK yielded similar PAHs content (3551 and 3412 ng mg -1 , respectively). Much higher PAHs content was observed for the Jet A-1/2,5-DMF blend (4662 ng mg -1 of soot). If we apply fuel sooting tendencies on the values presented in Table 4, we can obtain the following trend on total PAHs content for the tested fuel blends: Jet A-1/2,5-DMF > Jet A-1 > SPK > Jet A-1/DEC > Jet A-1/MOC > Jet A-1/1-Butanol. Therefore, excluding the Jet A-1/2,5-DMF blend, the addition of oxygenated biofuels to Jet A-1 reduces considerably the emissions of soot absorbed PAHs. The lowest PAHs content was obtained for Jet A-1/1-Butanol blend. Based on our results, the contribution of polyaromatics with different molecular weight (between 125-280 g mol -1 ) differs from one soot sample to another. Figure 6 gives the relative abundance of PAHs on soot; the size of the disks is proportional to the quantity of soot produced. Figure 6 shows that light PAHs (with 2 aromatic rings) represents only 3% or less of the total mass of PAHs quantified, while medium size PAHs (with 3-4 aromatic rings) constitute more than 58%. The highest concentration of heavy polyaromatics (≥ 5 aromatic rings) was observed for SPK (lowest sooting tendency) with 42% of the mass of soot collected whereas it reached only 28% of soot collected mass for the Jet A-1/2,5-DMF blend which has the highest sooting tendency. Global toxicity of soot samples was calculated based on the work of Nisbet and Lagoy [START_REF] Nisbet | Toxic equivalency factors (TEFS) for polycyclic aromatic-hydrocarbons (PAHs)[END_REF] that consider the individual toxicity of PAHs. The present work indicates (Table 5) that global toxicity of soot increases when Jet A-1 is blended with 1,2-DMF whereas for other blends, it decreases. Toxicity of soot particulates follows the trend: Jet A-1/2,5-DMF > Jet A-1 > SPK > Jet A-1/MOC > Jet A-1/DEC > Jet A-1/1-Butanol. Such impact of fuel blending, based solely on PAHs content, needs to be confirmed by direct test such as in vitro experiments.

CONCLUSION

In this work, we investigated the effect of blending conventional fuels (Diesel and Jet A-1) with oxygenated or synthetic fuels that can be derived from biomass. Carbonyl compounds emission and soot surface-bounded PAHs contents were measured using a four-cylinder Diesel engine and a premixed flame burner, respectively. The results indicated no Jet A-1/DEC correlation between oxygen content of the fuel blends and sooting tendency, PAHs concentrations, or carbonyls emissions. The global concentration of PAHs on soot was considerably lowered using oxygenated fuels, except 2,5-DMF. The lowest PAHs content was found for the 1-butanol-Jet A-1 blend. Conversely, the total carbonyl emission increased when oxygenated biofuels were added to the reference Diesel fuel. Ethanol and 1-butanol were found to increase considerably the emissions of carbonyl compounds. In all samples, acetaldehyde and pyrene were the most abundant compounds among carbonyl and polycyclic aromatic compounds, respectively. This study shows that blending conventional fuels with synthetic fuels does not imply a reduction of pollutants emission nor soot toxicity reduction.

ACKNOWLEDGMENTS

This

Figure 1 .

 1 Figure 1. Schematic representation of the CBCs sampling system used from engine measurements.

  Fuel O/C in fuel 1: Jet A-1 2: SPK 3: Jet A-1/1-Butanol 4: Jet A-1/DEC 5: Jet A-1/MOC 6: Jet A-1/2,5

Figure 2 .

 2 Figure 2. Schematic representation of the apparatus used for soot production and deposition.The deposited soot was removed mechanically, weighed, and introduced into 2 mL of HPLC-grade nhexane (Sigma Aldrich with stated purity ≥ 97%). The soot extracts were filtered twice using hydrophobic PTFE filter (pore size 0.2 µm by Alltech) to remove the non-soluble fraction mostly containing carbonaceous material. The resulting soot extracts were analyzed using a Shimadzu HPLC system leading to identification and quantification of PAHs and their derivatives[START_REF] Andrade-Eiroa | An alternative to trial and error methodology in solid phase extraction: an original automated solid phase extraction procedure for analysing PAHs and PAHderivatives in soot[END_REF]. The HPLC system was equipped with a C18 column (Grace Vydac 201TP, 250 mm × 4.6 mm × 5 µm), an automatic injector, and a photodiode array detector. Acetonitrile and water were used as mobile phase with the elution program: 0-30 min at 50:50 v/v followed by a linear gradient to 100/0 v/v at 90 min, and remained at 100/0 v/v until 110 min. The flow rate and column temperature were 0.5 mL/min and 30 °C, respectively. A certified 18 EPA PAHs standard (2000 µg/mL of each, Supelco) was used for identifying and quantifying naphthalene, 1-methylnaphthalene, 2-methylnaphthalene, acenaphthylene, acenaphthene, fluorene, phenanthrene, anthracene, fluoranthene, pyrene, benzo[a]anthracene, chrysene, benzo[b]fluoranthene, benzo[k]fluoranthene, benzo[a]pyrene, dibenzo[a,h] anthracene, benzo [g,h,i]perylene, and indeno[1,2,3-c,d] pyrene.

Figure 3 .

 3 Figure 3. Total emission of carbonyl compounds.

Figure 4 .

 4 Figure 4. Contribution of each carbonyl compound to total emissions of the fuels considered here.

Figure 5 .

 5 Figure 5. Individual carbonyl compounds emissions from the different fuels tested.

Figure 6 .

 6 Figure 6. PAHs relative abundance in the analyzed soot extracts. Small PAHs: MW = 128.17-166.22 g/mol; Medium PAHs: MW = 178.23-228.29 g/mol; Large PAHs: MW = 252.31-276.33 g/mol.

Table 1 .

 1 Fuel characteristics for engine test.

	Fuel	DCN	FSN
	DF DF/Ethanol * DF/1-Butanol * DF/DEC * DF/MOC * DF/TPGME * DF/SPK * DF/RME *	55.34 46.32 49.71 48.96 54.22 54.75 52.74 55.56	2.35 2.01 1.64 1.15 2.07 1.01 2.45 1.12
	* mixture of Diesel fuel/additive 90:10 v/v

Table 2 .

 2 Engine characteristics

	Number of cylinders	4
	Bore (mm)	76
	Stroke (mm) Displacement (cm 3 )	80.5 1460.74
	Number of injectors	4
	Compression Ratio	15.21

Table 3 .

 3 Fuel characteristics for burner experiments.

Table 4 .

 4 Summary of PAHs measurements.

	Fuel	Total 18 EPA PAHs (ng mg -1 of soot)
	1 Jet A-1	3551
	2: SPK	3412
	3: Jet A-1/1-Butanol	253
	4: Jet A-1/DEC	1978
	5: Jet A-1/MOC	819
	6: Jet A-1/2,5-DMF	4662

Table 5 .

 5 Global toxicity of soot samples*

	Fuel	Equivalent toxicity (TEQ)*	Variation to Jet A-1 %
	Jet A-1/1-Butanol	1.29	-99
	Jet A-1/DEC	10.83	-94
	Jet A-1/MOC	83.98	-57
	SPK	115.90	-40
	Jet A-1	193.57	0
	Jet A-1/2,5-DMF	574.14	+197
	* calculated global equivalent toxicity according to
	[30].

  work received funding from the European Research Council under the European Community's 7th Framework Programme (FP7/2007-2013) / ERC grant agreement n° 291049-2G-CSafe. The authors are grateful to J. Lemaire and F. Peyroux for technical assistance.

	NOMENCLATURE
	φ:	Equivalence ratio
	BMEP : Brake Mean Effective Pressure
	CBs:	Carbonyl compounds
	DCN: Derived cetane number, ASTM D7668
	DF:	Diesel fuel
	DNPH: 2,4-dinitro phenylhydrazine
	HPLC: High pressure liquid chromatography
	MW:	molecular weight
	PAHs: Polyaromatic hydrocarbons
	PTFE: Polytetrafluoroethylene
	RPM: rotation per minute
	SPK: Synthetic paraffinic kerosene
	UV:	Ultra-violet


L Jeannin

F Bignonnet
email: francois.bignonnet@univ-nantes.fr

F Agostini
email: franck.agostini@ec-lille.fr

Y Wang
email: yi.wang@ec-lille.fr

Stress effects on relative permeabilities of tight sandstones

Keywords: relative permeabilities, loading effect, sandstones, homogenisation perméabilités relatives, effet du chargement, grès, homogénéisation

This note deals with the study of stress-sensitive relative permeability experimentally observed in low permeable sandstones. These sandstones are made up of quartz grains surrounded by permeable interfaces between grains and pores. A micromechanical model of relative permeability behaviour under loading highlights the role of the closure of interfaces. Morphological models adapted to permeable sandstones show conversely that the relative permeabilities do not depend on the loading. To cite this article: A. Name1, A. Name2, C. R. Mecanique xxx (20xx)

Introduction

Low permeable gas sandstone reservoirs, also called tight reservoirs, are generally considered as stress-sensitive reservoirs [START_REF] Shanley | Factors controlling prolific gas production from low-permeability sandstone reservoirs: Implications for resource assessment, prospect development, and risk analysis[END_REF]. Numerous lab tests have shown that the absolute permeability of these reservoir rocks decreases strongly with confinement. For single phase flow, this dependence to confine-ment is related to the existence of joints and interfaces in tight rocks, which will close when loading is increased [START_REF] Fu | Experimental study of the stress dependence of the absolute and relative permeabilities of some tight gas sandstones[END_REF][START_REF] Ghabezloo | A micromechanical model for the effective compressibility of sandstones[END_REF][START_REF] Holditch | Tight gas sands[END_REF][START_REF] Schmitt | Characteriza-tion of Brazilian tight gas sandstones relating permeability and Angstrom-to micron-scale pore structures[END_REF][START_REF] Shanley | Factors controlling prolific gas production from low-permeability sandstone reservoirs: Implications for resource assessment, prospect development, and risk analysis[END_REF][START_REF] Walsh | The effect of pressure on porosity and the transport properties of rock[END_REF]. For two phase flow in porous media, not only the absolute permeability, but also the relative permeabilities can a priori be modified under loading. Recent experimental results demonstrate that, for tight sandstones, relative permeability curves will change with confinement [START_REF] Wang | Pétrophysique et micromécanique des grès tight en relation avec leur microstructure[END_REF].

The purpose of this paper is to use micromechanics in order to provide support to these experiments and to explain -at least qualitatively -why the behaviour of relative permeabilities under loading can be different between low permeable rocks and permeable rocks. The idea is to compare different morphological models representative either of permeable reservoir rocks or of tight rocks. This note is organised as follows. Section 2 presents a specific model describing tight rocks behaviour under loading. The proposed modelling is shown to be in fair agreement with measured data. Contrarily, considering permeable medium in section 3, we show for some simplified microstructures that the relative permeabilities should not depend on the loading.

Tight sandstones relative permeabilities under loading

Recent experimental results on low permeable sandstones [START_REF] Wang | Pétrophysique et micromécanique des grès tight en relation avec leur microstructure[END_REF] highlight a dependence of the gas relative permeability on the confinement (Figure 1). The tight sandstone studied by Wang et al. [START_REF] Wang | Pétrophysique et micromécanique des grès tight en relation avec leur microstructure[END_REF] is made up of quartz grains and pores (Figure 2); the contact zone between grains, partially cemented or not, forms joints with very thin openings e (i.e smaller than a few of micrometers) as compared to the grain radius R, which is around 100 µm (i.e. e R). From a fluid flow point of view, pores and interfaces define the pore network.

Grain with interfaces: the hydraulic model in the dry case

The morphological model chosen to represent sandstones is in agreement with the SEM observation concerning microstructure (Figure 2): tight sandstones are considered as an assemblage of spherical grains surrounded by permeable interfaces and of pores. The spherical shape is used to translate the fact that the considered phase does not manifest any anisotropy [START_REF] Dormieux | Homogenized models of stress-sensitive reservoir rocks[END_REF][START_REF] Kröner | Self-consistent scheme and graded disorder in polycrystal elasticity[END_REF].

We define the microscopic scale as the one of the heterogeneity of the microstructure (that is the scale of the grains with interfaces and of the pores). At the macroscopic scale, the representative elementary volume (rev) denoted Ω is regarded as a homogenised medium. The pore space occupies the domain Ω p in the rev. For tight sandstones, the porosity φ = |Ω p |/|Ω| is usually low (below 10%).

Let z denote the position vector at the microscopic scale in the rev. If a(z) is a field defined on Ω at the microscopic scale, its volume average is denoted by a = 1

|Ω| Ω a(z) dV . We also introduce the intrinsic averages a α = 1

|Ωα| Ωα a(z) dV over each phase α in Ω α , with α = p for pores or g for grains. Following their definitions, a as well a α are macroscopic quantities.

In order to derive an estimate of the macroscopic absolute permeability k hom and gas k rg and water k rw relative permeabilities, we will use a selfconsistent homogenisation scheme to relate microscopic fluid flow properties (i.e. at the pore scale) to the fluid flow behaviour of the rev at the macroscopic scale (see Appendix A).

As proposed for grains with permeable interfaces in [START_REF] Dormieux | Homogenized models of stress-sensitive reservoir rocks[END_REF], a two dimensional mathematical representation of the interfaces will be adopted. Note that we consider a viscosity equal to 1 in the following. At the microscopic scale, the core of a grain G i is assumed impermeable, while the surface velocity vector q is assumed to be proportional to the fluid pressure gradient in the interface:

z ∈ ∂G i : q(z) = -η grad |s p div |s q + [v] • n = 0 (1)
Following its definition, q lies within the tangent plane of ∂G i . The physical dimension of q (respectively η) is velocity times length (respectively permeability times length). Based on Poiseuille's law, the conductivity of the (half-)interface associated with one given grain is given by [START_REF] Dormieux | Homogenized models of stress-sensitive reservoir rocks[END_REF]:

η = 1 2 e 3 12 ( 2
)
where e is the opening of the interface.

In order to estimate the homogenised permeability k hom , two auxiliary problems are considered: one related to the flow in the pores and the other to the flow in the interfaces, around grains.

The first auxiliary problem is a classical Eshelby problem featuring a spherical pore embedded in an infinite homogeneous medium of permeability k hom (which is an unknown at this stage). At infinity, the boundary conditions on the fluid pressure write p(z) = ∇P 0 • z where ∇P 0 is an auxiliary macroscopic pressure gradient, which will be related later to the actual macroscopic pressure gradient ∇P . Since the characteristic size of the pores is large compared to the width of interfaces, the pore space is regarded as a phase with infinite permeability. The average of the pressure gradient and the velocity in the pore phase p are deduced from eqs. (A.1) and (A.2):

∇p p = 0 ; v p = -3k hom ∇P 0 (3)
In view of estimating the flow in the interfaces around the grains, we next consider a second auxiliary problem of a spherical grain G of radius R surrounded by its half-interface, embedded in an infinite homogeneous medium of permeability equal to the homogenised permeability k hom , with the same boundary conditions as previously. The resolution of this auxiliary problems leads to the following estimates of the pressure gradient and velocity averages over grains surrounded by an interface of thickness e:

∇p g(e) = 3/2 k hom k hom + η(e)/R ∇P 0 v g(e) = -3 k hom η/R k hom + η(e)/R
∇P 0 (4) using v g(e) = 3 4πR 3 ∂G qdS . Note that this result can be retrieved from eqs. (A.1) and (A.2) for a classical Eshelby problem by replacing the composite inclusion (grain+interface) by a uniform equivalent inclusion of permeability:

k eq (e) = 2η(e) R = e 3 12R (5)
Since all the interfaces in the material do not have the same thickness, we further assume a continuous distribution of interfaces by means of p.d.f. α; α(e)de represents the proportion of impermeable grains that have an interface with a thickness e in the range [e; e + de]. α is normalised such that:

∞ e=0 α(e)de = 1 (6)
From ( 4), ( 5) and the distribution α, the pressure gradient and velocity averages over all grains surrounded by interfaces is:

∇p g = ∞ e=0 3
k hom 2k hom + k eq (e)
α(e)de∇P 0

v g = ∞ e=0 -3 k hom k eq (e) 2k hom + k eq (e)
α(e)de∇P 0 (7) We now come back to the definition of the homogenised permeability k hom as the tensor relating the average over the rev of the pressure gradient ∇p to the average over the rev of the velocity v by v = -k hom • ∇p. Elimination of ∇P 0 in the estimates of ∇p and v via ( 3) and ( 7) yields the following estimate of the homogenised permeability [START_REF] Dormieux | Homogenized models of stress-sensitive reservoir rocks[END_REF]:

k hom = 3φk hom + (1 -φ) ∞ e=0 3 k hom k eq (e)
2k hom +k eq (e) α(e)de (1 -φ) e max /e min = 1 + where characterises the relative width of the distribution of the opening of grain joints, the solution to ( 8) is given by [START_REF] Bignonnet | Experimental Measurements and Multi-Scale Modeling of the Relative Gas Permeability of a Caprock[END_REF]:

k hom log-uni. = k eq (e med ) ρ 3 2 -ρ 1-3φ 1-φ -3 2 2(ρ 1-3φ 1-φ -1) (9)
If all interfaces have nearly the same width ( 1),

k hom log-uni. = k eq (e med ) 1-3φ 1 + 1+3φ 8(1-φ) 2 + o( 2 ).

Relative permeability model

Let us now evaluate the relative permeabilities in a simplified way. We denote the water saturation S w . To evaluate the gas relative permeability k rg (S w ), the previous hydraulic model is modified by assuming that interfaces and pores are now saturated either by water or by gas. Gas flow is assumed to take place only in pores and interfaces filled with gas, so that pores and interfaces saturated with water are impermeable with respect to gas flow.

In this system, the water phase constitutes the wetting fluid. We define e as the largest opening of the interfaces invaded by water: for a given saturation S w , each inclusion with interface opening e between 0 and e (S w ) is supposed to be filled up by water.

To sum up, the considered rev is made up of four different components:

-pores saturated with gas, with an infinite gas permeability, of volume fraction φ(1 -S w ), -pores saturated with water, impermeable to gas, of volume fraction φS w , -grains with an interface saturated with gas, of volume fraction (1 -φ) e=∞ e=e α(e)de -grains with an interface saturated with water, impermeable to gas, of volume fraction (1 -φ) e=e e=0 α(e)de The homogenised permeability is estimated as previously using the self-consistent scheme. By direct use of (A.6), the estimate of the homogenised gas permeability is the positive root to :

-φS w 1 2

pores w. water

-(1 -φ) e e=0 1 2
α(e)de grains w. water interfaces

+ φ(1 -S w )
pores w. gas

+ (1 -φ) ∞ e=e k eq (e) -k hom (S w ) 2k hom (S w ) + k eq (e)
α(e)de grains w. gas interfaces = 0 (10) when it exists, or zero. The transition between these two cases is obtained for the critical water saturation S w,c defined such that k hom (S w,c ) = 0 in [START_REF] Fu | Experimental study of the stress dependence of the absolute and relative permeabilities of some tight gas sandstones[END_REF], hence verifying:

(1 -φ) ∞ e=e (Sw,c) α(e)de + φ(1 -S w,c ) = 1 3 (11)
The relative permeability of the water phase k rw could be evaluated in a similar way. The end points of relative permeability curves take place in a natural way as the percolation threshold of the selfconsistent scheme, which may be not as accurate as dedicated percolation models as discussed later.

Unfortunately, we are not able to evaluate α in a straightforward way. The distribution β(e) of pores supplied by interfaces of largest thickness in the range [e; e + de] can be related to the capillary pressure curve 1 . The usual interpretation of capillary curve assumes that the filling of a pore is determined by a single interface thickness; in drainage, the pore associated to the largest interface gets filled by the non-wetting fluid at first and at increasing capillary pressures, pores fed with increasingly smaller interfaces become invaded by the non-wetting fluid. A simplified hypothesis consists in assuming that α(e) is proportional to β(e). Both distribution being normalised to the unity, we hence assume:

α = β (12)
e can be interpreted as a threshold of opening, below which all interfaces and connected pores are saturated in water. Under the assumption eq.( 12), e e=0 α(e)de = S w and the relative gas permeability is non-zero only for S w < 2/3 from eq. ( 11). Similarly, the relative water permeability is non-zero only for S w > 1/3. Gas and water phases may flow simultaneously for S w in [1/3; 2/3]. Note that the end points will be different for α = β or for a different morphological model.

Under the assumption [START_REF] Guéguen | Microstructures, percolation thresholds, and rock physical properties[END_REF] and in the case of a log-uniform distribution of interfaces, the solution to [START_REF] Fu | Experimental study of the stress dependence of the absolute and relative permeabilities of some tight gas sandstones[END_REF] and the definition k hom (S w ) = k rg k hom (0) lead to the following gas relative permeability estimate [START_REF] Bignonnet | Experimental Measurements and Multi-Scale Modeling of the Relative Gas Permeability of a Caprock[END_REF]:

k rg,log-uni. = ρ 3 2 -ρ 1-3φ(1-Sw ) 1-φ + 3 2 (2Sw-1) ρ 1-3φ(1-Sw ) 1-φ -1 ρ 1-3φ 1-φ -1 ρ 3 2 -ρ 1-3φ 1-φ -3 2 (13)

Loading effect

The purpose of this section is to present a poromechanical model of a low permeable sandstone consistent with the morphology of the previous section, in order to explain the relative permeability evolution under loading which has been observed experimentally. This constitutes an extension to relative permeabilities of a previous work in the dry case [START_REF] Dormieux | Homogenized models of stress-sensitive reservoir rocks[END_REF].

We hence turn to the micro-poromechanical model proposed by [START_REF] He | A poroelastic model for the effective behavior of granular materials with interface effect[END_REF]. The solid grains are assumed to be rigid and surrounded by elastic inter-faces. The constitutive law of the interface is characterised by a normal stiffness K n and a tangential stiffness K t , supposed to be the same for all interfaces. More precisely, for a macroscopic stress increment ∆Σ, let us denote ∆T the stress vector increment acting on the interface and [∆u] the displacement jump increment across the interface. The constitutive law of the interface reads:

∆T • n = K n [∆u • n] -∆p w or g ∆T t = K t [∆u t ] ( 14
)
where n denotes the normal vector to the interface,

∆T t = ∆T -(∆T •n)n and ∆u t = ∆u-(∆u•n)n.
Let us consider an isotropic stress increment ∆Σ = ∆Σ m 1. In what follows, we assume for simplicity that water and gas pressures do not vary significantly compared to the imposed confinement increment at the microscopic scale. Using the strain concentration rules established in the microporomechanical model of [START_REF] He | A poroelastic model for the effective behavior of granular materials with interface effect[END_REF], the variation of opening of each interface is proved to be independent of its initial opening e and given by:

∆e = K n R λ |∆Σ m | (15)
where λ is a morphological parameter which depends on the porosity. From a hydraulic point of view, a confining stress increment has thus a significant consequence as each interface with an initial opening lower than ∆e is closed and a part of the pore network is no more accessible to fluid flow. Since by definition β(e)de is the fraction of pores supplied by an interface of largest thickness in the range [e; e + de], the porosity ∆φ (≥ 0) trapped due to a load increment verifies2 :

∆φ = φ ∆e e=0 β(e)de (16)
To sum up, under the load increment the accessible porosity is φ-∆φ and the fraction of trapped pores ∆φ. Under assumption [START_REF] Guéguen | Microstructures, percolation thresholds, and rock physical properties[END_REF], the volume fraction of grains with open interfaces is then (1 -φ)(1 -∆φ φ ) and the volume fraction of grains with closed interfaces (1 -φ) ∆φ φ . The homogenised permeability is then estimated as previously using the self-consistent scheme. Due to loading, trapped pores or grains surrounded by closed interfaces are taken into account as being impermeable to any fluid flow. Let us denote by e (≥ ∆e) the opening of interface (before loading) which separates water filled interfaces from gas filled interfaces (after loading). By direct use of (A.6), the estimate of the homogenised gas permeability is the positive root to :

The function f (S w ) can be interpreted as a law of distribution of water in interfaces. Clearly, the function f has to increase from f (0) = 0 to f (1) = 1. In what follows, we will work under the simplifying assumption that the proportion of open, water filled interfaces among open interfaces is proportional to the water saturation, that is:

f (S w ) = S w (19)
Let us first analyse the effect of loading on the relative permeability end points. Under the assumptions ( 12) and ( 19), generalisation of (11) in order to account for closed interfaces and trapped pores based on [START_REF] Schmitt | Characteriza-tion of Brazilian tight gas sandstones relating permeability and Angstrom-to micron-scale pore structures[END_REF] indicates that the relative gas permeability vanishes for the critical water saturation:

S w,c = ∆φ -2 3 φ ∆φ -φ (20)
Hence, the end point of gas relative permeability under loading arises for a lower water saturation than for the unloaded case. As loading increases, the connectivity of the porous media decreases, and the critical water saturation for which gas relative permeability vanishes decreases. If the load closes more than two thirds of the interfaces, the permeability vanishes even in the dry case. In the vicinity of the critical saturation [START_REF] Walsh | The effect of pressure on porosity and the transport properties of rock[END_REF], the proposed model allows to retrieve the experimental observation that the gas relative permeability curve under loading is below the unloaded case one.

Further, under the assumptions ( 12) and ( 19) and in the case of a log-uniform distribution of interfaces, the homogenised permeability of the partially saturated medium under a load which traps a porosity ∆φ is found to verify:

k hom (S w ) = k eq (e med ) ρ 3/2 -ρ 3/2-2η+ζ 2(ρ η+ζ -1) (21)
where

η = 1 -∆φ φ (1 -f (S w )) ζ = ∆φ + (φ -∆φ)S w + (1 -φ) ∆φ φ (1 -φ) + (1 -φ) 1 -∆φ φ f (S w ) -2(φ -∆φ)(1 -S w )
(1 -φ) (22) under the approximation k eq (e -∆e) ≈ k eq (e) in the last integral of [START_REF] Schmitt | Characteriza-tion of Brazilian tight gas sandstones relating permeability and Angstrom-to micron-scale pore structures[END_REF] (which turns out to be a good approximation).

To study the influence of the shape of the p.d.f. α which governs the grain interface openings, a lognormal distribution of the opening with a mean of log(e) = -7 and a standard deviation of one has also been considered. The solution to eq. ( 17) is no longer analytical and has been computed numerically. The outputs of the model are presented in Figure 3 for both the log-normal and log-uniform distributions with a standard deviation of 1, in the case of an initial porosity of 0.1 and a trapped porosity ∆φ of 0.03. It is observed that our simple model does account for the stress-sensitivity of relative permeability curves for tight sandstones. At equivalent mean and standard deviation, the influence of the shape -log-uniform or log-normal -of the p.d.f. α of the size distribution of interface openings, is of second order on the relative permeabilities.

Limitations of the model

The main interest of the present model lies in the use of the same morphological model for the fluid flow and mechanical problems. This allows to explain qualitatively the effect of loading on relative permeabilities, at least close to the critical water saturation. However, let us stress out that self-consistent estimates do not have the accuracy of percolation models [START_REF] Guéguen | Microstructures, percolation thresholds, and rock physical properties[END_REF], realised for example on a 3D image of the pore network. In particular, the percolation threshold values arising from the selfconsistent scheme may not be accurate.

The present model further relies on a number of simplifying hypotheses, which allow us to reach an almost complete analytical treatment of the loaded problem. Namely, this simplified model does not account for complex multi-phase flow mechanisms such as snap-off or flow by film [START_REF] Dullien | Porous media, fluid transport and pore structure[END_REF]. Moreover, eq. ( 12) constitutes a strong assumption on the connectivity of the porous network. A closer comparison of our simplified model to the experimental results presented in Figure 1 shows some discrepancies for water saturation values which are not close to the critical water saturation (particularly in the case of log-uniform p.d.f. α): the model predicts that the relative gas permeability is higher in the unloaded case than in the loaded one, whereas experiments show the opposite behaviour.

Effect of the interface-bulk porosity distribution laws

The purpose of this section is to study the influence of the assumption [START_REF] Guéguen | Microstructures, percolation thresholds, and rock physical properties[END_REF] on the connectivity of the porous network. A relationship between α and β more general than ( 12) is now assumed:

where g is a function increasing from g(0) = 0 to g(1) = 1 since both distribution are normalised to the unity. The assumption (12) constitutes the particular case g(x) = x. Under the new assumption (23), the variation ∆e of interface opening due to loading is related to the variation ∆φ of porosity by:

∆e e=0 α(e)de = g ∆φ φ (24)
Recalling that for a given water saturation S w , f (S w ) denotes the proportion of water filled interfaces among interfaces still open after loading (see [START_REF] Shanley | Factors controlling prolific gas production from low-permeability sandstone reservoirs: Implications for resource assessment, prospect development, and risk analysis[END_REF]), the critical interface opening e which separates gas filled from water filled interfaces now verifies:

e e=∆e α(e)de = 1 -g ∆φ φ f (S w ) (25)
In the case of a log-uniform distribution of interfaces, the homogenised permeability of the partially saturated medium under a load which traps a porosity ∆φ under assumptions [START_REF] Shanley | Factors controlling prolific gas production from low-permeability sandstone reservoirs: Implications for resource assessment, prospect development, and risk analysis[END_REF] and ( 23) is analogous to [START_REF] Wang | Pétrophysique et micromécanique des grès tight en relation avec leur microstructure[END_REF], up to the minor modifications:

η = 1 -g ∆φ φ (1 -f (S w )) ζ = ∆φ + (φ -∆φ)S w + (1 -φ)g ∆φ φ (1 -φ) + (1 -φ) 1 -g ∆φ φ f (S w ) -2(φ -∆φ)(1 -S w ) (1 -φ) (26)
The effect of the choice of the laws f and g are illustrated in Figure 4 where the former assumption

f (x) = g(x) = x is compared to another (arbitrary chosen) assumption f (x) = g(x) = √
x. The latter choice, although arbitrary, corresponds to a case where the interfaces are water filled faster than bulk pores in imbibition. In that case, the sharp decrease with S w of the relative gas permeability observed for loaded tight sandstones is well reproduced.

Thus an appropriate choice of f and g should be motivated in future works by experimental observations or more sophisticated modelling approaches such as pore network models [START_REF] Bakke | 3-D Pore-Scale Modeling of Sandstones and Flow Simulations in the Pore Networks[END_REF][START_REF] Blunt | Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow[END_REF]. Nevertheless, at this stage, the model makes it possible to understand the influence of the closure of joints on the relative permeability in the case of tight sandstones.

f(x)=g(x)= √ x, unloaded f(x)=g(x)= √ x, loaded f(x)=g(x)=x, unloaded f(x)=g(x)=x, loaded

Model of relative permeabilities of permeable sandstones

Is is also interesting to use the previous approach with porous medium models adapted to conventional, permeable sandstones (i.e. non-tight). We assume pores modelled as spherical pores or tubes and show that in these cases relative permeabilities are at first order independent of the loading. For such permeable sandstones, the conventional approximation used in reservoir engineering and hydrogeology, which consists in considering the influence of the loading on the absolute permeability only and not on the relative permeabilities, thus seems reasonable.

Porous network made up of spherical pores

Unloaded case. We adopt the same approach as in section 2 and consider a porous medium medium made up of spherical pores and an impermeable solid phase. Following Markov et al. [START_REF] Markov | Permeability of the Fluid-Filled Inclusions in Porous Media[END_REF], a spherical pore of radius R embedded in a porous medium is equivalent to a spherical permeable inclusion of equivalent permeability3 :

k pore (R) = R 2 6 ( 27
)
Let us consider now as previously a statistical distribution of pore radius characterised by a probability α(R) that a pore has its radius in the range [R; R + dR]. From eq. (A.6), the homogenised absolute permeability verifies the following implicit equation:

- 1 -φ 2 + φ ∞ R=0 k pore (R) -k hom 2k hom + k pore (R) α(R)dR = 0 (28)
where k hom is non-zero only if φ ≥ 1 3 . This low percolation threshold is related to our choice of morphological model. Note that this feature could be avoided by choosing an equivalent cell permeability instead of spherical pores [START_REF] Boutin | Study of permeability by periodic and selfconsistent homogenisation[END_REF][START_REF] Bignonnet | Caractérisation expérimentale et micromécanique de la perméabilité et de la résistance des roches argileuses[END_REF] in which the pore volume is distributed around each grain.

Let us now turn to the estimation of the relative gas permeability. All pores whose radius is smaller than a critical size R governed by the saturation degree are assumed filled by water and impermeable to gas (S w = R R=0 α(R)dR). Using the selfconsistent scheme (A.6), the homogenised gas permeability k hom (S w ) = k hom (0) × k rg (S w ) verifies the following relationship:

φ ∞ R=R k pore (R) -k hom (S w ) 2k hom (S w ) + k pore (R) α(R)dR - 1 -φ + φS w 2 = 0 (29)
For this morphological model, the porosity φ must be greater than 1/3 and the critical water saturation for which the relative gas permeability vanishes is S w,c = 1 -1 3φ .

Loaded case. The variation of pore radius ∆R related to an isotropic mechanical loading increment ∆Σ m can be readily estimated from microporomechanics on the same morphological model (see [7] for more details). The resolution of the micro-poromechanical problem leads to:

∆R R = 1 3k sc (1 -α sc ) ∆Σ m (30)
where k sc and α sc are functions of the solid phase elastic properties and of the porosity [7]. It is further assumed that the pores do not close nor collapse due to loading, since obviously spherical pores are more difficult to close than to joints between grains. From ( 27) and (28), the variation of the gas permeability in the dry case is:

∆k hom (0) k hom (0) = 2 ∆R R = ∆k pore k pore (31)
However, equations ( 29) and (31) indicate that the relative permeability does not depend on the loading at first order, i.e. when changes of volume fraction of pores and grains related to loading are neglected.

Porous network made up of tubes

Unloaded case. Our second morphological model of conventional, permeable sandstones consists in a porous medium with an isotropic distribution of pores that are tubes (circular cylinders of radius r).

Let n be the unit vector along the axis of a tube. The permeability k n of a tube of radius r in its axial direction is given by Poiseuille's law:

k n (r) = γr 2 with γ = 1/8, (32)
while the permeability k t in directions normal to n is chosen equal to k hom [START_REF] Barthelemy | Effective permeability of media with a dense network of long micro fractures[END_REF]. The permeability tensor of such a tube is then the transverse isotropic tensor k n n ⊗ n + k t (1 -n ⊗ n). Using eqs. (A.1) and (A.2), the average of the pressure gradient and velocity in tubes of radius r whose axis is along the direction n is estimated by:

∇p p(n,r) = 1 • ∇P 0 v p(n,r) = -[k n (r)n ⊗ n + k t (1 -n ⊗ n)] • ∇P 0 (33
) Next, we consider an isotropic spatial distribution for the tubes and a probability α(r) that a tube has its radius in the range [r; r + dr]. From (A.6), the following self-consistent estimate for k hom is finally obtained:

k hom = 2φ 9 -7φ ∞ r=0 k n (r)α(r)dr1 (34)
The partially water saturated case is dealt with assuming that tubes with a radius r < r (S w ) are saturated with water and impermeable to gas. In that case, the gas permeability verifies:

k hom (S w ) = 2φ 9 -7φ + 8φS w ∞ r=r k n (r)α(r)dr1 (35) For this morphological model, the critical water saturation for which the relative gas permeability vanishes is S w,c = 1.

Loaded case. In order to deal with the effect of loading on the relative gas permeability, the variation ∆r of the radius of the tubes under an isotropic loading increment ∆Σ m can be assessed from micro-poromechanics on the same morphological model than for the permeability. From the hydraulic point of view, the relevant output of such poromechanical model is that ∆r/r is proportional to the load increment, i.e.: ∆r r = h∆Σ m (36

)
where h depends on the porosity as well as on the bulk and shear modulii of the solid phase. Combining ( 35) and (36), one concludes as for the previous model with spherical pores that the relative permeability is not sensitive to loading at first order.

Conclusion

The application of homogenisation techniques to simplified -but representative -morphologies of either tight, low-permeability sandstones or conventional, permeable sandstones allows to qualitatively retrieve an experimentally observed behaviour: the relative gas permeabilities of tight sandstones are much more sensitive to confinement than the ones of conventional sandstones. The underlying phenomena at the microscopic scale is the occurrence of closure of the joints due to loading in tight sandstones, which defines new percolation pathways. Conversely, it is shown that the relative permeability is not very sensitive to loading for conventional sandstones under the assumption that the load only modifies the pore radii, without modification of the connectivity of the pore network.

It is emphasised that the tight sandstone model has allowed us to identify that a required, key information is the connectivity of the interfaces and the bulk pores, through the functions f and g. The appropriate choice of these functions has not been investigated in the present work. It should be motivated in future works by experimental observations or dedicated modelling, for example using pore network models. Future works should also take into account capillary effects in the micro-poromechanical model. The difference in the pressure of the gas filled and water filled pores may indeed play an important role on the variation of joint openings, and thus affect the pore connectivity.

Figure 1 :

 1 Figure 1: Gas relative permeability measurements for 2 low permeable sandstone samples -loading of 3MPa (blue curve) -40 MPa (red curve)

Figure 2 :

 2 Figure 2: SEM image of a tight sandstone

∞ e=0 3 k 1 e

 31 hom 2k hom +k eq (e) α(e)de (8) In the case of a log-uniform distribution of interfaces α log-uni. (e) = ln(emax/emin) , if e ∈ [e min , e max ], denoting e med = √ e max e min and ρ =

  grains w. open, water filled interfaces+ (φ -∆φ)(1 -S w )open pores w. gas+ (1 -φ) ∞ e=e k eq (e -∆e) -k hom (S w ) 2k hom (S w ) + k eq (e -∆e)α(e)de grains w. open, gas filled interfaces = 0 (17) when it exists, or zero. To go forth, our model requires a key information: for a given water saturation S w , what is the proportion of water filled interfaces among interfaces still open after loading ? Let us denote by f (S w ) this fraction which reads, by definition of e : e e=∆e α(e)de = f (S w )

Figure 3 :

 3 Figure 3: Estimates of gas relative permeabilities without loading (line) and with loading (dashes) for φ = 0.1 and ∆φ = 0.03

Figure 4 :

 4 Figure 4: Estimates of gas relative permeabilities : effect of the laws f and g, for a log-uniform p.d.f. α with a standard deviation of 1, φ = 0.1 and ∆φ = 0.03.

which is measured without confinement of the rock sample.

Changes of volume fraction of pores and grains related to loading are here neglected.

Rigorously (27) is valid only in the limiting case where k pore (R) is large compared to the permeability of the surrounding porous medium.

Appendix A. Self-consistent estimate of the homogenised permeability

The homogenised permeability k hom is defined as the tensor relating the average over the rev of the pressure gradient ∇p to the average over the rev of the velocity v by v = -k hom • ∇p for appropriate boundary conditions (e.g. p(z) = ∇P • z on the boundary of the rev where ∇P = ∇p is the macroscopic pressure gradient).

The self-consistent homogenisation scheme used in this work relies on the estimate of the velocity and pressure gradient averages of each phase of the heterogeneous rev. These averages are estimated from the solution to Eshelby problems, devised as follows. For each phase, an inclusion I (or elementary particle representative of the phase) is embedded in an infinite medium whose permeability k 0 is uniform. An auxiliary uniform pressure gradient ∇P 0 is applied such that p(z) = ∇P 0 • z at infinity. If the inclusion I is an ellipsoid with homogeneous permeability k I , the pressure gradient and the velocity solution to the Eshelby problem are constant inside the inclusion. In this case, the pressure gradient inside the inclusion is given by:

• ∇P 0 (z ∈ I), (A.1) where P I 0 is the so-called Hill tensor of the inclusion I in the infinite medium of permeability k 0 . In the case where k 0 is isotropic (i.e. k 0 = k 0 1), the Hill tensors of a spherical and circular cylindrical inclusions are given by :

where n is the unit vector collinear with the axis of the cylinder. The velocity inside the inclusion is then simply deduced from the combination of Darcy's law with the Eshelby localisation rule (A.1) as:

3) The average of the pressure gradient and of the velocity over the rev are then estimated as :

The self-consistent estimate consists in choosing the permeability k 0 of the embedding infinite medium for each auxiliary Eshelby problem as the looked-for homogenised permeability k hom . The elimination of the auxiliary pressure gradient in (A.4) then provides the self-consistent estimate of the homogenised permeability as the solution to the equation:

where P hom (z) is the Hill tensor of the inclusion I comprising z embedded in the infinite medium of permeability k hom . The implicit equation (A.5) simplifies to:

where δk(z) = k(z) -k hom .