Florence Sagnard
email: florence.sagnard@ifsttar.fr

Jean-Philippe Tarel

Template-matching based detection of hyperbolas in groundpenetrating radargrams for buried utilities

Keywords: Ground-penetrating radar, pipe detection, hyperfrequency, template matching, soil characterization, ultra-wide band

Ground-Penetrating Radar (GPR) is a mature geophysical technique that is used to map utility pipelines buried within 1.5 m of the ground surface in the urban landscape. In this work, the template-matching algorithm has been originally applied to the detection and localization of pipe signatures in two perpendicular antenna polarizations. The processing of a GPR radargram is based on four main steps. The first step consists in defining a template usually from finite-difference time-domain simulations made of the nearby area of the hyperbola apex associated with a mean size object to be detected in the soil whose mean permittivity has been previously experimentally estimated. The second step consists in a pre-processing on the raw radargram to correct variations due to antennas coupling, then the template matching algorithm is used to detect and localize individual hyperbola signatures in an environment containing unwanted reflections, noise and overlapping signatures. The distance between the shifted template and a local zone in the radargram based on the L1 norm allows to obtain a map of distances. A user-defined threshold allows to select a reduced number of zones having a high similarity measure. In a third step, in each zone minimum or maximum discrete amplitudes belonging to a selected hyperbola curve are semi-automatically extracted. In the fourth step, the discrete hyperbola data (i,j) are fitted by a parametric hyperbola modeling 1 1 2 3 4 5 6 7 8 9 1 2

using a non linear least square criterion. The algorithm has been implemented and evaluated on numerical radargrams, and afterwards on experimental radargrams.

Introduction

Ground-penetrating radar (GPR) is a well-known non destructive technique for imaging shallow subsurfaces by a propagating electromagnetic energy downward into the ground from a transmitting antenna that is reflected at subsurface boundaries between media possessing different electromagnetic (EM) properties and is collected by a receiving antenna. GPR is a proposed technology to map utility pipelines in urban environments (fiber optics, telecommunication lines, electrical cables, water and gas pipes, district heating network), most of them have been buried within 1.5 m of the ground surface [START_REF] Benedetto | Civil engineering applications of ground penetrating radar[END_REF][START_REF]Proceedings of the General Meeting of COST Action TU1208[END_REF]. The buried targets produce characteristic hyperbolic signatures in the radargram (Bscan) issued from the moving of the GPR system along a linear path. Hyperbolic signatures represented in the distance-time domain are specific to the target size, shape and dielectric characteristics and the orientation of the electric field [START_REF] Radzevicius | Ground penetrating radar polarization and scattering from cylinders[END_REF][START_REF] Huici | Accurate ground penetrating radar numerical modeling for automatic detection and recognition of antipersonnel landmines[END_REF]. Challenges lie in detecting and classifying targets in an environment with variations in the surface cover (asphalt, paving, sand, gravel, grass), in the subsoil with spatial vertical or horizontal variability of the soil texture (water content, backfill soil, multilayers…), and the presence of buried targets that can be close to each other (overlapping hyperbolas) and have wide variations in dimensions, and in dielectric properties (metallic or non-metallic).

An ultra-wide band (UWB) ground-coupled radar operating in the frequency band [0.46 ; 4] GHz and made of bowtie slot antennas has been preferred to an air-launched radar because it increases energy transfer of electromagnetic radiation in the sub-surface and penetration depth [START_REF] Benedetto | Civil engineering applications of ground penetrating radar[END_REF][START_REF] Sagnard | Design of a compact ultra-wide band bow-tie slot antenna system for the evaluation of structural changes in civil engineering works[END_REF]. Moreover, UWB GPR and SFCW (step frequency continuous wave) have been used to probe the soil structure using the benefit from both low and high frequencies bringing a compromise in terms of depth resolution and penetration in a single measurement. This paper is focused on the processing and analysis of radargrams using the semi-automatic templatematching algorithm applied to GPR application in order to recognize hyperbolas produced by buried targets and extract information from them: detection, localization, and characterization. Special attention has been paid to study the effects of the dielectric contrasts between the targets and the surrounding soil and also of the polarization configuration of the GPR system relative to the main axe of canonical targets (pipes or strips).

In the literature, several works deal with the semi-automatic buried target detection and characterization. There are two main classes of methods: in the first class the model of the pattern is given a priori or designed by hand, in the second class the model is learned from the collected data. The methods based on Hough transform and its derivatives [START_REF] Windsor | A data pair-labeled generalized Hough transform for radar location of buried objects[END_REF][START_REF] Wang | Fast detection of GPR objects with cross correlation and Hough transform[END_REF][START_REF] Huang | A Automatic detection of buried utilities and solid objects with GPR using neural networks and pattern recognition[END_REF], genetic algorithm (GA) [START_REF] Kobashigawa | Classification of buried targets using ground penetrating radar: Comparison between genetic programming and neural network IEEE Antennas Wireless Propag[END_REF], and image segmentation and hyperbola fitting [START_REF] Shihab | Radius estimation for cylindrical objects detected by ground penetrating radar[END_REF] are in the first class.

The methods based on artificial neural networks (ANNs) [START_REF] Huang | A Automatic detection of buried utilities and solid objects with GPR using neural networks and pattern recognition[END_REF][START_REF] Gamba | Neural detection of pipe signatures in ground penetrating radar images[END_REF][START_REF] Kobashigawa | Classification of buried targets using ground penetrating radar: Comparison between genetic programming and neural network IEEE Antennas Wireless Propag[END_REF][START_REF] Nunez-Nieto | GPR signal characterization for automated landmine and UXO detection based on machine learning techniques Remote[END_REF], support vector machine (SVM) [START_REF] Muniappan | Utility detection and parameter estimation using ground penetrating radar[END_REF], and wavelet analysis [START_REF] Cieszczyk | The curvelet transform application to the analysis of data received from GPR technique Electrika IR[END_REF][START_REF] Zhan | GPR measurement of the diameter of steel bars in concrete specimens based on the stationary wavelet transform[END_REF][START_REF] Ni | Buried pipe detection by ground penetrating radar using the discrete wavelet transform Computers and[END_REF] are in the second class. The template matching [START_REF] Jain A K, Duin | Statistical pattern recognition: A review[END_REF][START_REF] Vernon | Machine vision-Automated visual inspection and robot vision[END_REF][START_REF] Zyada | 3D template based automatic landmine detection from GPR data[END_REF][START_REF] Jain A K, Duin | Statistical pattern recognition: A review[END_REF] is widely used in pattern recognition, since it is a simple method, based on a prior model extraction from one or several images, which leads to a good accuracy for detecting targets with reduced variability in observed images. Surprisingly, this approach has not been previously used in GPR, when observed patterns due to objects of interest are of quite similar aspects in the radargrams.

The template-matching algorithm is readily implemented on a computer and has a reasonable calculation time providing that the image and the template sizes are not too large. The template does not need to contain a hyperbola response very similar to the one to be detected, it only needs to be able to discriminate between hyperbolas of interest and the background (noise and small heterogeneities in the soil compared to the targets). The analysis of the distance map based on the L1 norm using an optimal threshold has allowed to select a limited number of template positions. In this work, the use of two antenna orientations along a profile (named polarization diversity) has allowed to obtain two distance maps from which a mean distance map should improve hyperbola detection in one of the two radargrams where the polarized response appears too weak and thus not detectable. The fitting of each hyperbola applied on maximum or minimum amplitudes is performed by a parametric analytical model based on the straight ray path hypothesis and using the constraint least square (LS) criterion.

Considering canonical objects (pipes and strips), the template-matching algorithm has been tested on simulated (commercial software Empire) and experimental radargrams. This paper is structured as follows: in section 2, the radar geometry is briefly described that was modeled using FDTD simulations. The parametric analytical ray path modeling was used to model hyperbolas. Section 3 deals with the steps of the pattern recognition algorithm that includes mainly the template-matching method and the parametric hyperbola fitting. The developments and validation were first performed on synthetic images. In section 4, a few experimental images acquired in a laboratory site are analyzed and the results associated with the parameter analysis are discussed. Finally, in section 5 conclusions and perspectives are drawn.

Modeling of the GPR system

A. Numerical FDTD modeling

The GPR SFCW system made of a pair of shielded bowtie slot antennas designed on a FR4 substrate, with real relative permittivity ε '=4. 4 and thickness e=1.6 mm (see Fig. 1a), has been preliminary studied using FDTD simulations (commercial software Empire) [START_REF] Sagnard | Design of a compact ultra-wide band bow-tie slot antenna system for the evaluation of structural changes in civil engineering works[END_REF] on the soil surface on a distance close to 1.6 m appears prohibitive for a parametric study (more than one week using a CPU i7-950 3.07 GHz), a pair of more simple planar antennas, named blade dipoles working in the same frequency band, but with a higher reflection coefficient, has been used. The blade dipoles designed on a FR4 substrate are non shielded and their dimensions are 290×56 mm 2 (see Fig. 1c) [START_REF] Tebchrany | Assessement of statisticalbased clutter reduction techniques on ground-coupled GPR data for the detection of buried objects and cracks in soil[END_REF][START_REF] Rejiba | Zero-offset profiling using frequency cross-hole radars in a layered embankment test-site: antenna design[END_REF].

The offset between antennas in simulations and experiments has been fixed to 60 mm, and the elevation h s above the soil is equal to 10 mm. The soil electrical parameters ( ε s ' , σ s ) are assumed constant across the frequency range. The GPR system is moved linearly on the soil surface with a step Δy=40 mm (see Fig. 1a) to acquire a radargram. Because the experimental and modeled antennas appear as complementary UWB dipoles, the polarization of the electric field ( E θ , E ϕ ) appears perpendicular. In the simulations, the excitation current has the shape of the derivative of the Gaussian function with a time zero estimated to 0.3 ns and a duration (99% of the total energy) of 0.5 ns (peak frequency 1 GHz, bandwidth 3 GHz).

B. Analytical ray path modeling

The analytical modeling based on the ray path hypothesis supposes that the target is localized in the far-field zone of the GPR system. Such a modeling helps the interpretation of the earliest reflection component of a hyperbola pattern. Considering a buried canonical object, such as a cylindrical pipe, the linear displacement of the GPR system on the soil surface (h=0) such as depicted in Fig. 1a gives a time-distance curve with the shape of a hyperbola. The two-way (round-trip) travel time can be expressed as a function of the horizontal position y of the radar, the radius R , its horizontal location at y 0 and its depth d under the soil surface [START_REF] Shihab | Radius estimation for cylindrical objects detected by ground penetrating radar[END_REF]. The lateral center-to-center antenna distance SR is considered in this study because of the marked dimension size of the antennas in the two main polarization configurations.

The equations associated with the travel-time write as follows:

y T = y i -SR / 2 y R = y i + SR / 2 ¿ { ¿ ¿ ¿ ¿
(1)

T Tx 2 t arg et = [( y 0 -y T ) 2 +( d + R ) 2 ] 0 . 5 -R T t arg et 2 Rx = [( y 0 -y R ) 2 +( d + R ) 2 ] 0 . 5 -R ¿ { ¿ ¿ ¿ ¿
The velocity v of the medium depends on the real relative permittivity ε s ' such as:

} } ≪ε rSub { size 8{s} } rSup { size 8{'} } \) } { ¿ v=c/√ε s ' ( ε s ¿ ¿ (2)
The generalized hyperbola equation ( y i , t i ) including the antenna offset is expressed by:

t i = ( T Tx 2 t arg et +T t arg et 2 Rx ) / v (3)
Thus, the hyperbola depends on five parameters ( SR , y 0 ,d, R ,v ).

A preliminary parameter study presented in Fig. 2 was conducted considering two configurations of the SFCW GPR system: parallel (SR=291 mm), and end-fire (SR=422 mm) [START_REF] Sagnard | Design of a compact ultra-wide band bow-tie slot antenna system for the evaluation of structural changes in civil engineering works[END_REF]. Assuming a soil with ε s ' =3. 5 ( s =0) and a buried pipe at depth d=180 mm with several radius values [0,40,80,120,160] mm, we first remark from Fig. 2a that when SR=0 the radius reduces the arrival time as if the velocity was increased; when R=0, a time difference of 0.64 ns at the apex is observed between the configurations SR=0 and SR=291 mm. We observe that the antenna distance SR flattens the apex and produces a delay in the arrival time, particularly significant when SR=422 mm (see Fig. 2b). The time difference at the apex increases with the soil permittivity as visualized in Fig. 2c. In the configuration SR=291 mm, we remark that the arrival time difference at the apex with a pipe radius less than 40 mm appears less than 0.

Template matching technique and hyperbola fitting

The proposed algorithm based on template-matching aims to detect semi-automatically hyperbola signatures without a preliminary training period, and to estimate target characteristics as reliably as possible.

The radargram processing is made of a series of steps:

-A pre-processing performed on the raw radargram to remove the antenna direct coupling and the clutter removal to enhance radargram quality;

-The construction of an amplitude distance map based on the translation on the radargram of a predefined template at every possible positions where a mean amplitude distance is evaluated according to the L 1 norm [START_REF] Huang | A Automatic detection of buried utilities and solid objects with GPR using neural networks and pattern recognition[END_REF][START_REF] Gamba | Neural detection of pipe signatures in ground penetrating radar images[END_REF]. A threshold value allows to select local discrete minima on the distance map that corresponds to a high level of similarly with the template.

-At the selected positions, the hyperbola points close to the apex that can correspond to a maximum or a minimum amplitude are extracted. For each curve, a fitting of the points is performed using a parametric model;

-The estimation of the parameters describing the hyperbola curve such as d ,R , v is obtained using a hyperbola fitting of the points to the analytical relation (3) according to the LS criterion. The Hessian matrix in 2D and particularly its eigenvalues can be used to characterize the uncertainties on the estimated parameters (d ,v ) .

A. Pre-processing

In a radargram, the strongest and first signal detected associated with antenna groundcoupling may mask the hyperbolas of shallow buried objects. Depending on the depth of the objects and more precisely the degree of interaction of their responses with the signal of the clutter, an adequate clutter reduction technique has to be chosen [START_REF] Benedetto | Civil engineering applications of ground penetrating radar[END_REF][START_REF] Rejiba | Zero-offset profiling using frequency cross-hole radars in a layered embankment test-site: antenna design[END_REF]. In this study, we have considered objects buried at a depth higher than 10 cm, thus the median subtraction at each time sample is assumed sufficient to reduce the horizontal component of the clutter.

Preliminary to image processing, the duration of the radargrams considered has to be limited (10 ns for numerical data, and 5 ns for experimental data) in order to eliminate the soil background without potential targets. Moreover, the distance-time units ( Δy=40 mm in numerical and experimental data) of the template image defined have been scaled according to the image to be analyzed that supposes to perform a time interpolation; as in synthetic radargrams the time step was fixed to Δt=2.396 10 -11 s and 1.791 10 -11 s in the parallel and end-fire configurations, respectively (see Fig. 3a and3b), in experiments it was fixed to 5.5610 -11 s . A time zero correction has been applied in each radargram according to the shape of the excitation signal (the first positive amplitude component arrives here at 0.5 ns) to obtain a time zero at the wave ground bounce.

It must be underlined that in a radargram the data representing discrete amplitudes can be either binary or signed, depending on the quality and the dynamic range before the conversion; as in simulated data there is no background noise or unwanted reflections, the conversion of signed amplitudes into binary amplitudes doesn't change the image quality. In experimental data, the signed amplitudes will be preferred.

B. Template-matching technique

The template-matching algorithm supposes to preliminary define a template including the average aspect of the different hyperbolas of interest in the radargrams to be analyzed. A preliminary study involving buried pipes under study was made to visualize and interpret a few synthetic hyperbola signatures. This study is an extension of the work in [START_REF] Radzevicius | Ground penetrating radar polarization and scattering from cylinders[END_REF][START_REF] Sagnard | Using polarization diversity in the detection of small discontinuities by an ultra-wide band ground penetrating radar[END_REF]. Synthetic radargrams have been obtained from FDTD simulations (software Empire) considering a system with blade dipoles (offset=60 mm, h s =10 mm) in both configurations, parallel and endfire (see Fig. 3a and3b) moving linearly on a homogeneous soil ( ε s ' =3. 5 , σ s =0.01 S .m -1 , h=1 cm ). A pipe (radius R=32 mm) with fundamental dielectric characteristics such as air- filled, conductive, and dielectric ( ε ' =15 ) was buried at depth d=168 mm [START_REF] Radzevicius | Ground penetrating radar polarization and scattering from cylinders[END_REF][START_REF] Huici | Accurate ground penetrating radar numerical modeling for automatic detection and recognition of antipersonnel landmines[END_REF]; the depth d is measured between the soil surface and the top surface of a pipe. The hyperbola signatures visualized in Fig. 4a, 4b, 4c, 4d, 5a 5b and 5c are plotted in the time range 8 ns, and the time vectors were interpolated to match the experimental time step Δt=5.56 10 -11 s .

In general, we remark that the parallel configuration such as presented in 4a and 4b (SR=116 mm) generates multiple reflections with higher amplitudes in the hyperbola legs whatever the dielectric characteristics of the pipe as compared to the end-fire case. These reflections correspond to wave bounces between the soil surface and the pipe top. In the case of a non conductive pipe, waves propagate in the pipe, and we observe a marked reflection at the bottom of the pipe when it is filled with a medium having a high contrast with the soil, in particular in water, that induces a wave velocity difference. Moreover, the amplitude of the backscattering waves appears reduced when the soil is attenuating.

In the end-fire configuration, such as visualized in Fig. 4c and 4d (SR=340 mm), hyperbola patterns have few components, and appear more compact. The conductive pipe induces a specific inverse polarity (negative amplitude) in the reflection pattern that allows to distinguish it. In both polarizations, we remark that the conductive pipe has the strongest reflection amplitude. In the case of a high dielectric pipe ( ε ' =15 ), Fig. 5a and 5b highlight that in both polarizations multiple reflections induced between the top and bottom of the pipe can be distinguished because the velocity is relatively low inside the pipe; the parallel polarization gives a higher amplitude response. In Fig. 5b, artifacts that are induced by the clutter removal treatment appear visible because in the present case the clutter component and the target signature overlap (see Fig. 5c). Samples of these synthetic pipe signatures can serve as template images in the experiments.

The template does not need to be perfect, it only needs to allow us to discriminate between hyperbolas of interest and the background (small heterogeneities). When the set of hyperbolas of interest is too large to be represented by a single pattern, several patterns can be used in sequence. Each defined template includes a small portion of a hyperbola pattern in the vicinity of the apex and is scaled in time, distance and amplitude to the radargram under analysis.

Assuming an observed image g(i, j) with a size M ×N , and a template image t(i , j) , both with scaled amplitudes, we define the L 1 norm distance map E between g and t by the following relation [START_REF] Vernon | Machine vision-Automated visual inspection and robot vision[END_REF]:

E(m , n)= ∑ i ∑ j |t (i, j )-g(i-m , j-n )| (4)
The summation is evaluated at all pixels (i,j) of the template t that is translated to all possible positions (m,n) along the observed image g. This leads to an image of L1 distances, also named absolute distance map. The position (m,n) at which the smallest value E(m ,n) is obtained corresponds to the best match between the template t and the corresponding subimage in g. A threshold value allows to select a limited number of local minima corresponding to distances E less than the threshold where the template is well matched with enough amplitude (visualized by "+" signs on the radargram, see Fig. 7c and8c).

The template-matching algorithm has been first applied on synthetic data. The GPR system made of blade dipoles in the parallel (see Fig. 6) and end-fire configurations has been considered on a homogeneous soil structure including three close buried pipes at several depths. The pipes, from the left to the right, made of conductor (n°1), dielectric (n°2, ε s ' =9

) and air (n°3, ε s ' =1 ) are separated by a distance equal to 300 and 250 mm respectively. The pipe radii are respectively 32, 22 and 11 mm. The radargrams after clutter removal are presented in Fig. 7c and 8c in the parallel and end-fire configurations (time range 10 ns). The comparison of both radargrams highlights a difference in shape and amplitude in the hyperbola signatures of the three objects that is explained by excitation direction of the electric field towards the pipes which are either conductive or dielectric. In both polarizations, we observe that the air-filled pipe has a significant lower response than the two other pipes, and thus appears masked; this phenomenon is caused by the weak permittivity contrast between the pipe and the soil that will be higher in a real soil.

In the end fire configuration, the synthetic template ( 189×7 pixels, Δt=1.79 1 10 -11 s ) visualized in Fig. 8a represents the upper part, near the apex, of the hyperbola signature of the conductive pipe in the end-fire configuration to further consider hyperbolas that overlap because of close objects; the apex appears in the middle of the time range. The complementary version of this template ( 141×7 pixels, Δt=2.396 10 -11 s ) has been used in the parallel configuration as it appears more adapted. At present, the template includes a significant upper zone without signal to allow the detection of first hyperbolas containing information.

In the parallel configuration, the similarity measure, applied on Fig. 7c (considering a time range 8 ns to reduce late reflections) while translating the template of Fig. 7a, leads to the distance map of Fig. 7b. In this last figure, two marked minima can be visualized that correspond to the positions of both significant hyperbolas. Defining the maximum threshold value to 120, the first two local positions of hyperbolas are detected, as visualized in Fig. 7c (time range 10 ns). We remark that the air-filled pipe has not been detected. In the end-fire configuration, the analysis of the radargram presented in Fig. 8c using the template of Fig. 8a according to the template matching algorithm, gives the distance map of Fig. 8b. The maximum threshold value of 162 allows to detect the two first hyperbola zones associated with the conductive and dielectric pipes such as visualized in Fig. 8c (time range 10 ns). As previously, the pipe air-filled pipe cannot be detected.

The polarization diversity (joint information from both perpendicular polarizations) can be used to enhance the detection in each polarization [START_REF] Sagnard | Using polarization diversity in the detection of small discontinuities by an ultra-wide band ground penetrating radar[END_REF]. Thus, a mean distance map can be calculated from the individual distance maps in the parallel and end-fire configurations. A validation with experimental data of this point is described in Section 4.

C. Hyperbola fitting

For each detected template location, the next processing step consists in three sub-steps: location of the maximum or minimum amplitude arch associated with first arrival times, extraction of the discrete amplitudes, and fitting of the arch by an analytical hyperbola equation.

The curve points associated with the first arrival times of a hyperbola may correspond to a maximum or a minimum amplitude. Because higher order reflections in a hyperbola pattern may produce a stronger amplitude as compared to the amplitude of the first reflection, an user interaction is necessary (semi-automatic) to select a hyperbola curve either on the upper or on the lower half zone of the template position. Starting from the middle point of the template, close points belonging to the hyperbola curve on the left and on the right legs (usually 3 or 4 points) are extracted step-by-step. It is possible to extract points outside the template location mainly when no signature overlaps. Because the number of curve points extracted is limited, a polynomial fitting of the second order was made to refine the estimate of the abscissa y 0 at the apex. Then, knowing a priori parameters SR and y 0 , a fitting with the analytical modeling described by equation ( 3) and based on a constraint LS criterion has allowed to estimate the three parameters (d ,R, v) . The constraint consists in defining for each parameter a minimum value, such as in the present case 40 mm (or 30 mm), 50 mm and 70 mm.ns -1 for R, d and v, respectively and a maximum value 60 mm, 200 mm and 190 mm.ns -1 , respectively.

In the analysis of the numerical radargrams, initial values are respectively 40 mm, 100 mm and 134.16 mm.ns -1 (for ε '=5 ). In order to confirm the stability of the algorithm, a second fitting in a similar manner based on the analytical formula is performed using the last estimated parameters of R, d and v as initial values. Each hyperbola fitting in the radargram is made independently of each other.

Considering the radargrams in Fig. 7c and 8c associated with the parallel and end fire configurations, the parameter values obtained from the hyperbola fitting are collected in Tables 1 and2, respectively. In each template location, the hyperbolas fitted were plotted (see Fig. 7c and8c). From the parameter values, we remark that the pipe radius cannot be properly estimated, as it appears smaller than the radiating aperture of an antenna (on the order of the first Fresnel zone in the soil [START_REF] Rial | Resolution of GPR bowtie antennas: An experimental approach[END_REF]); thus, the radius estimate remains equal to the initial value defined in the optimization function. In general, the real permittivity value was estimated with an uncertainty less than 15%, and this uncertainty influences the estimation of the depth. We observe that the curve fitting of pipe n°2 gives a lower permittivity estimate in both configurations as compared to the curve fitting of pipe n°1. In the end-fire configuration, we remark from Fig. 8c that the pipe signature n°2 is made of a significant second reflection that appears quite higher that the first reflection, that makes it difficult to detect. Consequently, the second hyperbola signature which arrives 0.4 ns later has been analyzed and the depth will need an adjustment of -69.3 mm (see Table 2). The fitting of this signature gives a higher objective function fval for pipe n°2 because 4 points has been considered in the left leg of the hyperbola. Otherwise, the objective functions are of the order of 10 In the sandy box, a 25 mm diameter dielectric pipe and a thin vertical 10 mm width (2 mm thick) conductive strip were buried at depths estimated to 160 and 170 mm respectively such as presented in Fig. 9. Both objects are separated by a 750 mm distance. Firstly, a synthetic template was computed from 3D FDTD simulations (software Empire) considering the detailed bowtie slot antenna geometries (section 2.A) and 32 mm radius conductive pipe buried in a soil with a real relative permittivity ε '=3 .5 ( σ=10

-2 S .m -1
). In this template visualized in Fig. 10, it was wise to define a compact hyperbola signature with reduced multiple reflections to further detect different hyperbolas in a radargram in both polarizations.

Thus, a conductive pipe in the end-fire configuration has been considered. The template (see The template-matching algorithm was performed on the experimental radargrams visualized in Fig. 11b and 11d and associated with both parallel and end-fire polarizations. The L1 norm distance maps are presented in Fig. 11a and11c. The maximum threshold values leading to the detection of the first and most significant hyperbolas are 0.162 and 0.195 respectively;

higher values give additional detections (false alarms) that don't correspond to buried objects but to background heterogeneities. We remark that the end-fire polarization does not permit to detect the air-filled pipe. The results of the parameter evaluation from the LS fitting of each hyperbola detected are collected in Tables 3 and4. In general, the positions y 0 of the objects appear correctly evaluated. Concerning the real permittivity value of the soil, the parallel configuration gives higher estimates as compared to the end-fire configuration, and consequently the target depths appear more important; from Table 3, we remark that the depths of the pipe and the strip have been both evaluated to 200 mm, and more important than those a priori evaluated (see Fig. 9). In the case of a soil having weak permittivity variations, an additional step would be to find an optimum permittivity value issued from the several estimates. We remark that the objective function fval associated with the LS fitting is presently higher for experimental data than for synthetic data, of the order of 10 -2 , because the image quality appears lower.

To detect the air-filled pipe in Fig. 11d, joint information from radargrams issued from the parallel and end-fire configurations (polarization diversity) could be used. Thus, the mean distance map was calculated from the individual distance maps in the parallel and end-fire configurations (see Fig. 11a and11c) leading to the result visualized in Fig. 12a. A weak threshold (0.18) situated between the previous ones (0.162 and 0.195) associated with both polarizations, has allowed here to detect both hyperbolas without false alarms in the end-fire configuration, that was not possible when this polarization was only considered.

Further insight into the solutions of the parameters estimated in the fitting has been gained by calculating the Hessian matrix H at the stationary point to evaluate its nature using its eigenvalues.

The rate of convergence and sensitivity to round-off errors is given by the condition number of matrix H, that is the ratio of its largest to its smallest eigenvalues. In the present examples, fixing the pipe radius that cannot be evaluated properly has led to a decrease of the condition number. The eigenvalues associated with the several fitting in both polarizations are collected in (4.0 10 -5 ; 9.6 10 -3 ) Table 5: Eigenvalues of the Hessian matrix at the estimates of the depth d and velocity v from hyperbola fitting in both polarizations A second soil structure built in the sandy box includes two buried objects, a 25 mm diameter water-filled pipe and a horizontal 2 mm width strip with a 10 mm height such as visualized in Fig. 13a. In the parallel configuration, the radargram not presented here shows two very similar weak hyperbola signatures quite difficult to detect. In the end-fire configuration, the radargram presented in Fig. 13c shows two different hyperbola signatures. Using the template of Fig. 10, the template-matching algorithm has permitted to localize only the strip signature (threshold value 0.202) because the water filled pipe shows a quite different and specific hyperbola pattern; a threshold value higher than 0.202 will lead to detect several signatures in the background. The parameter estimations issued from the hyperbola fitting lead to the values collected in Table 6

Conclusion

This study is focused on the feasibility of detecting the location, depth, lateral dimension of long but thin buried objects such as pipes or strips in a soil using a GPR system. The most difficult tasks in the analysis of GPR radargrams is the split of overlapping signatures, the detection of a hyperbola pattern in a noisy background and in a weak image quality, the detection of an object with a small lateral dimension and a weak contrast with a perturbed surrounding medium. The work presented was carried out in two main steps: firstly, a theoretical parameter study (based analytical and numerical results) has consisted in analyzing hyperbola signatures of different dielectric pipes as a function the dielectric contrast with the surrounding medium, the polarization of the electric field and their diameter, and secondly a semi-automatic algorithm based on the template-matching technique and LS hyperbola fitting was developed to detect, localize and characterize a given but not restrictive hyperbola pattern. The template-matching algorithm does not need a training period, but needs the help of the user because of the great diversity of hyperbola signatures and the objective of the algorithm is to operate with a non perfect template.

The template-algorithm validated on numerical radargrams issued from FDTD simulations has shown that distinct templates have to be associated with both polarizations. Depending of the level of details required, several templates may be successively used to identify 20 39 progressively the dielectric nature of a buried target among several in a radargram. A template has to consider the lateral dimension of the GPR system when it is quite larger than the size of the object to be detected. A sensitivity analysis has to be perform to study in details the influence of the soil permittivity used in the template to detect hyperbolas in a radargram. The advantage of using a SFCW GPR is that the excitation signal used in the simulated template can be used in the IFFT transform. Even with a few extracted hyperbola points, the present work shows that main parameters associated with the target can be estimated, but may be slightly correlated. This correlation can be diagnosed by means of the Hessian matrix.

Considering measurements in a controlled environment, this work has shown that synthetic templates can be used to analyze radargrams provided that the templates have the same time step and amplitude range adjusted if necessary by interpolation; thus, a set of patterns has to be collected. The analysis of an experimental radargram in the end-fire configuration has shown that in the case of a buried water-filled pipe and a vertical conductive strip, the selected synthetic template was not able to detect the signature of the water-filled pipe because of its specific signature. The template-matching algorithm applied to GPR was extended by using the benefit of polarization diversity. This implies to define a mean distance map using the distance maps of both polarizations to improve the detection of weak hyperbola signatures in a polarization and strengthen the algorithm robustness. Further studies will be focused on bringing improvements to the template-matching algorithm, particularly during its validation on experimental radargrams acquired in the laboratory site Sense-City [26] where different targets representing utilities and cracks were buried. Thus, we are planning to build a data base of templates with variable sizes and variable dielectric contrasts of pipes and strips and variable surrounding materials. This new site is a good opportunity to refine the clutter removal procedure and to evaluate the performances of the template matching algorithm in a situation close to a real environment.

LIST OF FIGURES

3 .

 3 . The antenna radiation characteristics are adequate for operating in an UWB in the frequency range [0.46; 4] GHz very close to the soil surface (ground-coupled). Each antenna is enclosed in a shielded conductive rectangular box filled with a three-layered absorbing foam with dimensions 362× 231 × 67.5 mm Because the simulation time of the pair of antennas moving

4 .

 4 Evaluation on experimental dataTargets with different dielectric characteristics have been buried in a large sandy box of the public square Perichaux, Paris 15 th district. The sand was not compacted and its relative dielectric permittivity was estimated to 3.5 from CMP (Common Mid Point) measurements.The time zero was corrected according to CMP measurements that have allowed to locate the direct waves. The pair of bowtie slot antennas (see Fig.1b) separated by a 60 mm offset was moved linearly on the soil surface with a step of 40 mm to acquire a radargram. Both polarizations, parallel and end-fire, have been considered separately.The complex transmission coefficient S21 at the receiving antenna port was measured and stored using a vector network analyzer (VNA) ANRITSU MS 2026B in the frequency range [0.05;4] GHz (1601 samples). To obtain the transmitted signal in the time domain, an apodization (zero padding) has been applied to S21 in order to smoothly extend the frequency band from 4 to 9 GHz. The excitation signal used in FDTD simulations, the first derivative of the Gaussian function with a duration equal to 0.5 ns, has been introduced in each experimental signal in the frequency domain. The product of the spectrum of the excitation signal (convolution in the time domain) with the measured S21 has been performed to further calculate the inverse Fourier transform (IFFT) and obtain time data.

Fig. 10 )

 10 Fig. 10) was scaled in time and amplitude to match the experimental time step

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 12345 Figure 1: (a) Scheme of the GPR system imaging a subsurface with a buried pipe; (b) Geometry of the bowtie slot antenna used in the measurements; (c) Geometry of the blade dipole used in the simulations

Figure 6 :Figure 7 :Figure 8 :

 678 Figure 6: Geometry of the GPR system, used in FDTD simulations, considering of a pair of blade dipoles moving on a soil including 3 different buried pipes (n°1 metal, n°2 dielectric

Figure 9 :

 9 Figure 9: Geometry of the sandy box with two buried objects, an air-filled pipe and a horizontal conductive strip, probed experimentally using a pair of bowtie slot antennas

Figure 10 :Figure 11 :Figure 12 :

 101112 Figure 10: Synthetic template issued from FDTD simulations involving a pair of bowtie slot

  Δt=5.56 E-11 s used in experiments. With a larger antenna distance SR=422 mm, the radius variations appear more easily separable at the apex, however the hyperbola apex appears more flat and its amplitude will appear lower when the arrival time is higher. The variation of the soil real relative permittivity in the range[4 ; 10] (R=40 mm) and illustrated in Fig 3chighlightsthat the hyperbola slope increases with the permittivity value, and the apex appears sharper.

	(at the apex when R=0, if SR=291 mm t0=2.89 ns, and if SR=422 mm t0=3.45 ns); thus, in
	this range the radius value cannot be distinguished in a radargram and particularly when
	considering a time step of
	1 ns

Table 1 :

 1 -3 . In general, we observe that the uncertainty of the depth is less than 22%. Parameter estimation in the parallel configuration for the synthetic Bscan of Fig.7c

		d	R	'	t 0	x 0	fval
	Pipe n°1	184.3 mm	40	3.25	2.3 ns	-189.1	1.68 10 -3
	(1 st
	maximum)					mm
	Pipe n°1	168 mm	32	3.5		-185 mm
	(true
	values)
	Pipe n°2	98.8 mm	40	2.92	1.26 ns	113.1 mm 7.67 10 -4
	(1 st
	maximum)
	Pipe n°2	79 mm	20 mm	3.5		115 mm
	(true
	values)
		d	R	'	t 0	x 0	fval
	Pipe n°1	150 mm	30	3.93	2.94 ns	-188.3	8.23 10 -3
	(1 st
	minimum)					mm
	Pipe n°2	149.9 mm	30	2.99	2.55 ns	113.9 mm 1.67 10 -2
	(2 nd
	minimum)	80.6 mm
	(1 st minimum
	by				2.11 ns
	extrapolation
	)

Table 2 :

 2 Parameter estimation in the end-fire configuration for the synthetic Bscan of Fig.8c

Table 5 .

 5 In general, the condition number is high, that means that correlation may exist between the two parameters and thus the convergence of the estimation algorithm appears

	slow.
		d	R	'	t 0	y 0	fval
	Pipe n°1	200 mm	60	3.6	3.08 ns	458.9	2.68 10 -2
	(maximum
	)					mm
	Pipe n°1	~160 mm	12.5 mm	3.5-4		~500 mm
	(true
	values)
	Strip n°2	200 mm	60	3.47	2.99 ns	1232 mm 2.23 10 -2
	(maximum
	)
	Strip n°2	~170 mm	5 mm	3.5-4		~1200
	(true
	values)					mm

Table 3 :

 3 Parameter estimation in the parallel configuration for the experimental Bscan of Fig.

				11b
		d	R	'	t 0	y 0	fval
	Strip n°2	165.6 mm	25.5	2.65	2.82 ns	1250 mm 1.85 10 -2
	(maximum
	)

Table 4 :

 4 Parameter estimation in the end-fire configuration for the experimental Bscan of Fig.

	11d

  . The permittivity estimate appears close to the value issued from the previous experiments in the end-fire configuration.

		d	R	'	t 0	y 0	fval
	Strip n°2	187.8 mm	60	2.49	2.61 ns	1299 mm 1.07 10 -1
	(maximum
	)

Table 6 :

 6 Parameter estimation in the end fire configuration for the experimental Bscan of Fig.

	13c

11


Muhammad Kashif
email: kashif091@gmail.com

Jérôme Bonnety
email: jerome.bonnety@upmc.fr

Alexis Matynia
email: alexis.matynia@upmc.fr

Patrick Da Costa
email: patrick.da_costa@upmc.fr

Guillaume Legros
email: guillaume.legros@dalembert.upmc.fr

Sooting propensities of some gasoline surrogate fuels: Combined effects of fuel blending and air vitiation

Keywords: gasoline surrogate fuels, yield sooting index, soot volume fraction maps

come

Introduction

Soot particles formed along combustion processes can be ultimately emitted into the atmosphere. Even at low concentration, the presence of these particles in urban air poses a serious public health problem [1,2] and may play a critical role in global warming due to a significant capacity to absorb the incoming solar radiation [3]. Therefore, strategies aiming at the reduction of soot emission by combustion devices are required.

While diesel engines exhibit some relatively high performances, especially in terms of CO 2 specific emission (g/km) [4], they suffer from the acknowledged dilemma called NO

x -soot trade-off. Indeed, soot and NO x emissions cannot be reduced at the same time by the diesel combustion process itself.

The technological solution currently implemented is to reduce NO

x formation rate inside the cylinders. To this end, using exhaust gas recirculation (EGR) lowers combustion temperature and inhibits nitrogen chemistry [5]. In contrast, soot particles that form in the rich regions of the fuel spray experience a reduced oxidation rate due to the lower temperature but they can be captured and burnt downstream in the exhaust flow through a particulate filter.

However, significant amounts of soot, especially ultra-fine particulate matters, can still escape from these soot traps [6]. Alternative technologies like Gasoline Direct Injection (GDI) have shown to increase combustion efficiencies exhibited by spark-ignited (SI) engines, therefore can help reduce CO 2 specific emissions by these engines [7]. However, this kind of strategies also leads to the formation of rich region in the vicinity of the fuel spray, which in turn increases soot formation, therefore emission from gasoline engines.

Due to these technological limitations, the strategies aiming at the reduction of soot emission by reciprocating engines need to be complemented by methodologies of fuel formulation [8,9]. The guidelines formulated by the recent Worldwide Fuel Charter [START_REF]Worldwide Fuel Charter, ACEA, the Alliance of Automobile Manufacturers, EMA, JAMA, and[END_REF] support these methodologies as the recommended aromatics content of gasoline has for instance been reduced from 50 % down to 35 % within few years.

Complex and competing physico-chemical processes are involved in soot formation during combustion [START_REF] Frenklach | [END_REF]. The comprehensive understanding of in-ception, growth and subsequent oxidation of soot particles remains a hot issue, especially with respect to complex blends of hydrocarbons [3,12]. More and more studies have focused on soot formation/oxidation for the last two decades, contributing to bridge the knowledge gap required to optimize fuel formulation [3,13,14]. Due to the recent development in computational power, refining the understanding of the aforementioned soot related processes can be achieved when contrasting experimental measurements with data extracted from numerical simulations of reference flames [5,15]. This method then requires sophisticated computational and experimental tools.

At atmospheric conditions, many practical fuels are liquid mixtures of hundreds of hydrocarbons which have variable and coupled effects on regulated engine emissions and especially soot production [3,9]. To develop a molecular-level understanding of soot formation/oxidation processes along combustion of these practical fuels, many works employ surrogate fuels and gradually increase the complexity of fuel formulation of these surrogates. Thus, special efforts have been devoted to the design of surrogate fuels for gasoline, diesel and jet fuels [16]. As an illustration of relevant surrogate fuels, blends of n-heptane, isooctane and toluene [16,17] have been shown to fairly reproduce the combustion behavior of gasoline [18]. Among other characteristics, the propensity of a surrogate fuel to produce soot then needs to be assessed [14,19]. To compare the sooting tendencies of a large number of chemical compounds, surrogate fuels, hydrocarbon based commercial fuels and mixture of hydrocarbons and bio-fuels, different methods like the smoke point height, the threshold sooting index (TSI), and the Yield sooting Index (YSI) are measured in laminar diffusion flames [20][21][22][23][24]. Among these protocols, the YSI methodology has lately been shown to be especially relevant and provide with fairly low uncertainties as the sooting propensity is then inferred directly from soot volume fraction measurements in the flame [22,24]. This sooting tendency is to be quantified in flames where the fuel to be investigated burns in homogeneous inflows. In this way, the discrimination among physical and chemical processes leading to soot formation in more practical devices will be allowed. The step that the determination of the fuel's YSI represents can then be considered a valuable contribution to the knowledge required to develop strategies for the reduction of the soot release by combustion devices.

Studying the sooting propensities exhibited by Primary Reference Fuels (PRFs), Kashif et al. [24] highlighted the blending effect that the evolution of YSI as a non-linear function of the isooctane mole fraction in the PRF reveals. Furthermore, these authors showed that a non-negligible cross effect of the isooctane mole fraction and CO 2 dilution of air on YSI exists. Thus, linear-by-volume (lbv) models to estimate the combustion properties, such as YSI, of these simple bi-component fuels would give poor predictions. This behaviour is attributed to synergistic and antagonistic effects whose magnitudes depend not only on fuel blending but also on air vitiation. These insights make the assessment of the combined effects of fuel blending and air vitiation crucial for surrogate fuels. Indeed, due to the relatively simple formulation of a surrogate fuel, the effect of fuel blending can be identified, therefore modeled. Moreover, investigating the effect of air vitiation by CO 2 dilution is of practical relevance as CO 2 is a major component of EGR.

In the present study, the experimental methodology recommended by

McEnally and Pfefferle is reproduced in order to measure the sooting tendencies of some gasoline surrogate fuels, i.e. binary mixtures of n-heptane/toluene and isooctane/toluene. Indeed, Yan et al. [25] noticed that for mixtures of n-octane and toluene, a parabolic curve fits the evolution of the experimental sooting tendencies with toluene mole fraction in the mixture (see Fig. 6 in Ref. 25). Later studies then showed that blending toluene and isooctane also leads to that kind of non-linear effect [17,19] while blends of toluene and n-heptane show fairly linear trends. In the following, the consistency of the YSI methodology is therefore assessed. The aforementioned blending effect is identified in terms of YSIs exhibited by n-heptane/toluene and isooctane/toluene mixtures. The methodology extended by Kashif et al. [24] to the sooting propensity in vitiated air is then applied. This enables the characterization of the combined effects of blending and CO 2 dilution of air. Eventually, as a contribution to the experimental database that the current computational efforts request, the whole set of measured soot volume fraction fields measured in the reference flames established is delivered as a supplemental material.

Quantification of sooting tendencies

Experimental methodology

The experimental setup used to reproduce the methodology prescribed by McEnally and Pfefferle [23] has been carefully described in Ref. 24. The diffusion flames are established over an axis-symmetric coflow burner identical to the one described by Santoro et al. [26] and used in previous investigations [24,27,28]. In the following, the axis of symmetry is (Oz) and its origin is located at the burner tip, defining the height above the burner (HAB). The cross-stream coordinate is r, which is the distance from the axis of symmetry.

In their original work [23], McEnally and Pfefferle suggested to evaluate the apparatus-independent Yield Sooting Index (YSI) of a fuel from the maximum soot volume fraction f v,max = max {f v (r, z meas ) | r 2 [0, R]} measured in the non-premixed flame burning in coflowing air. The fuel stream is composed of methane which is nitrogen diluted and doped at a given concentration with the vapor of the fuel whose sooting propensity is to be evaluated.

In the original methodology, the measurement f v,max is performed at a fixed height z meas above the burner. Then the fuel's YSI is calculated using Eq.

(1):

Y SI = C f v,max + D (1)
Here, C and D are apparatus-specific parameters that are determined by the arbitrary YSI values attributed to two reference fuels. For the present study, YSI values of benzene (30) and toluene (43.5 ± 3%) reported by

McEnally and Pfefferle [23] are used as references. YSI benzene is also the lower reference in the study of McEnally and Pfefferle [23]. YSI toluene is here chosen as the upper reference as it will be the sootiest fuel investigated along the present study. YSI toluene was delivered by McEnally and Pfefferle with an uncertainty of 3%. The mean YSI toluene is selected here which does not affect the interpretation of the results.

In the following, three similar methods provide with alternative evaluations of f v,max that the two-dimensional fields of soot volume fraction allow. As prescribed by Kashif et al. [24], implementing these different methods enables the assessment of the YSI consistency. Thus, Y SI Max. Ax. fv , Y SI Global M ax. f v and Y SI Max. Int. fv are the measurements of the sooting propensity when f v,max is the maximum soot volume fraction on the flame axis, the peak soot volume fraction into the whole distribution, and the maximum value of the integrated soot volume F v (z), respectively. The latter distribution is given by the following expression:

F v (z) = 1 ⇡R 2 Z R 0 2⇡rf v (r, z)dr (2)
where R=5.5mm is the inner radius of fuel tube.

The Light Extinction Measurement (LEM) technique as implemented by

Kashif et al. [27] provides with the two-dimensional soot volume fraction fields f v (r, z) required to evaluate the different kinds of f v,max and input them in Eq. (1). To this end, the ray that outcomes from a 100 mW Continuous

Wave laser operating at 645 nm is expanded to a collimated beam with a diameter of 70 mm. After the beam crossed the flame, it is decollimated. A pinhole with a diameter of 800 µm ± 5 µm is located at the focal point to provide a telecentric configuration possessing depth invariant magnification.

The pinhole also allows the slight beam steering induced by the temperature gradient to be mitigated as this steering can bias the deconvolution process.

The absorption of the beam through the flame is imaged by a Photon Focus MV1 12-bit progressive scan monochrome camera whose lens is equipped with a narrow band filter centered at 645 nm (±2 nm) and with a band width at one half the transmissivity maximum of 20 nm. A digital pulse generator controls the occurrence and the duration of the camera's CMOS sensor exposure, together with the opening of a shutter that chops the laser source. A frame grabber records the frames captured by the camera on a computer. As a line-of-sight technique, LEM is combined with a subsequent deconvolution that incorporates a Tikhonov regularization [27]. The local extinction coefficient distribution  ext (r, z) can then be computed from the projected data, i.e. the field of extinction coefficient integrated along the lineof-sight. This LEM technique has been shown to provide both fine temporal and spatial resolutions [27].

As prescribed by Zhou et al [29], the soot particles are here supposed to experience a relatively short residence time in the flame, therefore exhibit a range of diameter that is small enough to warrant the Rayleigh assumption. Thus, the local soot volume fraction is finally inferred from the local extinction coefficient as follows:

f v (r, z) =  ext (r, z) /(6 ⇡ (1 + ↵ sa ) E(m)) (3)
where E(m) is a function of the complex refractive index m of soot. ↵ sa is the relative contribution of scattering to extinction, therefore is directly related to soot morphology.

Thus, a careful choice of the soot refraction index is crucial to the consistency of the methodology. While the values of E(m) (about 0.28 at 645 nm) reported by Krishnan et al [30] are almost identical for n-heptane, isooctane, toluene and benzene, a significant discrepancy of ↵ sa appears. In the following, the values of ↵ sa recommenced by Krishnan et al [30] are selected, i.e. 0

for n-heptane and isooctane, and 0.1 for benzene and toluene. This lead to values of the denominator in Eq. ( 3) equal to 5.3 for n-heptane and isooctane, and 5.8 for benzene and toluene. For binary mixtures of n-heptane/toluene and isooctane/toluene, ↵ sa is supposed to follow a linear dependency on the mole fraction of toluene in the blend. This adjustment is required to take into account the increasing fraction of light scattering in the phenomenon of light extinction by the soot particles when toluene mole fraction is increased.

Indeed, with increasing toluene mole fraction, the particles exhibit higher characteristic diameters due to both hastened soot formation processes and longer residence time [30], as evidenced here (see section 3.1) by early soot load in the flame and increased flame length.

Nonetheless, ↵ sa stands in the denominator in Eq. ( 3). Accounting for the relative contribution of ↵ sa or not then induces discrepancies in YSI of the order of the experimental uncertainties. Finally, the aforementioned adjustment of ↵ sa hardly affects the calculated YSI values and the overall YSI trends.

Procedure

The procedure that is extensively documented in Ref.24 has also been followed along the present study. The reactants came from high-purity gas cylinders (CH 4 and N 2 ; 99.9 % stated purities), reagent-grade bottles (nheptane, isooctane and toluene; 99.5 %), and a compressor (air).

The only difference with the procedure reported in Ref. 24 is the composition of the carrier gas which is now a mixture of methane (50% in volume) and nitrogen (50% in volume), as opposed to pure methane used in the aforementioned study. The dilution of the carrier gas with nitrogen is here required to keep the sootier flames investigated below the smoke point conditions.

To evaluate the YSI of any blend stored in the tank of the setup, the carrier gas flow is first set through the Bronkhorst Controlled Evaporator and Mixer (CEM) and a lightly sooting diffusion flame of the carrier gas is established over the burner. Once the temperature of the CEM heater and the heated lines from the CEM to the burner is stabilized well above the condensation temperature of the liquid blend vapor, the liquid blend flowing from the tank through a Coriolis mass flow controller is vaporized. Its vapor is injected into the carrier gas. The flame then becomes brighter as soot production is clearly enhanced. A 10 s long sequence of LEM is recorded every minute at a frame rate of 30 Hz and with an exposure time of 15 ms, showing that any raw extinction field exhibits a fluctuation of its peak lower than 5% over 10 minutes after the liquid flow was established. After the last recording, the liquid flow rate is set to zero. The flame then progressively turns blue, showing that for any condition the fuel stream is only composed of the carrier gas, one minute after the liquid flow rate is set to zero. The carrier gas flow rate is also set to zero. Once the flame is extinguished, the liquid fuel tank is flushed and the whole line from the tank to the burner tip is purged with an inert gas, i.e. Argon, during 5 minutes to remove all the liquid and vapors from the system. The tank is then refilled with the following liquid blend to be investigated.

The experimental parameters that were kept constant along the present study are reported in Tab. 1, together with those set by McEnally and Pfefferle. The values 30 and 43.5 were here assigned to the YSIs of benzene and toluene, respectively. The temperature along the fuel line, i.e. 150 C, is similar as that set by McEnally and Pfefferle. The coflowing air flow rate, i.e. 60000 cm 3 /min, is higher here as this level was required to produce stable flames within the whole range of conditions investigated. However, some experiments conducted with a fuel line temperature of 180 C on one hand, and with a coflowing air flow rate of 80000 cm 3 /min on other hand, lead to YSI modifications lower than 0.5. Thus, the influence of these parameters is considered negligible as compared to the standard deviation among the YSIs provided by the different methods described earlier.

The parameters varied in the present study are:

• the composition of the bi-component blend, identified by the toluene mole fraction X T ol. in the mixture of toluene and n-heptane or toluene and isooctane. This parameter is adjusted when mixing different quantities of n-heptane or isooctane with toluene before refilling the tank.

• the vapor mole fraction X vap,i (i=1,2) into the fuel stream. This is controlled by the liquid flow rate through the Coriolis mass flow controller. The mass flow rate required to set a given X vap,i is calculated with ±1% uncertainty from compiled liquid phase densities and species molecular weights [31].

• the carbon dioxide mole fraction in the oxidizer stream, referred as X CO 2 . A mass flow controller allows the carbon dioxide flow rate to be adjusted. This flow is then mixed with the constant air flow rate to complement the oxidizer flow.

For every parameter, the range investigated within the frame of the present study is specified in Tab. 2, together with the corresponding estimated uncertainty.

The dopant concentrations investigated correspond to two reference benzene mass flow rates, i.e. 0.42 mg/sec (1.5 g/hr) and 0.56 mg/sec (2.0 g/hr), at room temperature. These flow rates correspond to the vapor mole fractions X vap,1 =1.86 10 2 and X vap,2 =2.47 10 2 , respectively. The mass flow rate for a given liquid fuel to be investigated is then adjusted such that its vapor mole fraction is equal to X vap,i (i=1,2), calculated with an absolute uncertainty of approximately 0.1% from compiled liquid phase densities and species molecular weights [31]. The reported uncertainties of the mass flow controller used for liquid mixtures are of the order of 0.1% of measured value. The uncertainties related to the use of compiled densities and molecular weight are conveyed into the YSI measurements as these quantities are required to compute the mass flow rates of the different mixtures.

In the following, any YSI reported is the mean value among the 150 evaluations allowed by the 150 fields of soot volume fraction measured for every set of conditions.

Effects of fuel blending and air vitiation

Both single effects, i.e. fuel blending on one hand, and air vitiation on the other hand, are first assessed. Combined effects are then addressed.

For every set of conditions investigated, the experimental database, i.e. every distribution of f v (r, z) in the whole flame, is delivered as a supplementary material. of the CMOS sensor, the detection threshold of the methodology in terms of peak soot volume fraction is here evaluated at 0.15 ppm [24].

Effect of fuel blending

For both kinds of binary mixtures and at any CO 2 content of the oxidizer, the peak soot volume fraction in the doped methane diffusion flame increases monotonically with increasing toluene mole fraction in the blend. Indeed, with increasing toluene mole fraction, the C/H ratio increases. Still, the soot formation rate is governed (at the first order) by the fuel C/H ratio [32].

Furthermore, the addition of aromatic compounds hastens soot formation in flames [19]. Therefore, increasing toluene mole fraction in the fuel stream results in a higher peak soot volume fraction. This also leads to a longer flame since the total fuel consumption is delayed due to the relatively slow soot chemistry.

To further elucidate the effect of blending on soot formation, the evolution of the integrated soot volume fraction F v (z) in the doped flames is plotted against the height above the burner in Figs. 2 for both (a) n-heptane/toluene and (b) isooctane/toluene mixtures. The first trend which also reveals qualitatively in Figs. 1 is the upstream location of the first soot particles detected which shifts gradually further upstream in the flame with increasing toluene mole fraction at any CO 2 content of the oxidizer. Concomitantly, the peak integrated soot volume fraction increases monotonically with the increasing amount of toluene into these blends while its location also shifts upstream.

Increasing toluene mole fraction then leads to higher slopes of the increasing parts of the curves which indicates an increasing soot formation rate. This results in larger soot volume fractions at any given HAB. These trends are qualitatively in agreement with the available data delivered by both coflow [14] and counterflow [19] burners.

However, a discrepancy between both mixtures appears in the increase rate of sooting propensity with toluene mole fraction. As revealed in Fig. 2(b), the peak soot volume fraction, therefore the soot formation rate, exhibited by the mixture of isooctane and toluene increase first rapidly with toluene mole fraction (0  X T ol.  0.4 at X CO 2 =0%). Then this increase starts decaying (X T ol.

0.6 at X CO 2 =0%). On the opposite, the increase exhibited by the mixture of n-heptane and toluene appears to be more uniform.

To further support these comments, YSIs are calculated following the three alternatives prescribed in section 2.1. Figure 3 exhibits the evolutions of YSI as functions of toluene mole fraction in the blends of n-heptane/toluene (blue dots) and isooctane/toluene (black dots) for the three CO 2 contents of the oxidizer investigated. The YSI scale is defined by reference values of 30 and 43.5 ± 3% which are attributed to benzene and toluene, respectively, for X CO 2 =0%. For a matter of legibility, all mean experimental YSIs, delivered by the three alternatives at both vapor mole fractions investigated, are not discriminated. These evaluations have been shown to be self-consistent as they lead to very similar evolutions of YSIs [24]. Second-order polynomials decently fit these experimental measurements. The procedure to derive them is given in section 3.3. The solid lines shown in Fig. 3 represent the evolutions of these polynomials. For both mixtures of n-heptane/toluene and isooctane/toluene, YSI increases monotonically with increasing toluene mole fraction, following the tendency exhibited by the absolute measurements shown in Figs. 1 and2. The aforementioned discrepancy in the increase rate of sooting propensity with toluene mole fraction also reveals in terms of YSI.

On one hand, an almost straight line can readily fit the calculated YSIs of the n-heptane/toluene mixture at X CO 2 =0% with reasonable accuracy showing a linear evolution of this mixture's sooting propensity as a function of toluene mole fraction. This linear trend turns into a slightly non-linear one at higher CO 2 contents of the oxidizer. Thus, the sooting propensity seems to be affected by a coupling between fuel blending and CO 2 dilution effects. On the other hand, a parabolic curve fits more properly the evolution of YSI exhibited by the isooctane/toluene mixture. This non-linear evolution of the sooting tendency with increasing concentration of toluene mole fraction can be explained in the light of the chemical mechanisms responsible for the formation of soot precursors species, i.e. poly-aromatic hydrocarbons (PAHs). Indeed, a synergistic effect of blending on the evolution of sooting propensity (measured here as YSI) has already been unveiled for binary mixtures of isooctane and toluene [19,25].

As shown in Fig. 33 in Ref.3, the pyrolysis of n-heptane produces mainly short chain alkenes -especially ethylene-while isooctane preferentially leads to longer chain alkenes (propene, butene, pentene, and hexene). These latter species then especially decompose to form propargyl radicals (C 3 H 3 ) and vinylacetylene (C 4 H 4 ). As an illustration, among the reactions leading to the formation of the aforementioned radicals, the decomposition of propene

(C 3 H 6 ) to propargyl radicals (C 3 H
3 ) via a series of H-atom abstractions is considered the preferential kinetic pathway [3].

In addition, Choi et al. [19]

C 6 H 5 CH 3 + H ⌦ A1 + CH 3 [3, 34, 35] (4)
A1 CH 2 + C 3 H 3 ⌦ A2 + H 2 [34] (5)
A1 C 6 H 5 + C 4 H 4 ⌦ A2 + H [35] (6)
Thus at lower toluene mole fraction the relatively high production of H radicals by heavy hydrocarbons such as n-heptane or isooctane may boost benzene and naphthalene formation, resulting in an increased sooting tendency that would have been under-predicted by a simple linear-by-volume model. Further enhancing this effect, the higher concentrations of propargyl C 3 H 3 radicals and vinylacetylene C 4 H 4 especially produced along the decomposition of isooctane hastens naphthalene formation and its growth to heavier PAHs.

The influence of these pathways among the overall soot formation processes is then reduced at higher toluene mole fraction due to the relatively lower concentrations of species originated from the pyrolysis of isooctane, resulting in relatively lower concentrations of propargyl radicals and vinylacetylene. As a consequence, the synergistic effect of blending on soot formation weakens with increasing mole fraction of toluene in isooctane/toluene mixture.

Effect of air vitiation

The effect of CO 2 dilution of the coflowing air stream reveals readily in Figs. 1. The trends match the observations already reported in the literature related to methane diffusion flames doped with vapors of liquid fuels established over the Santoro's burner [14,24,33]. The sooting propensity of doped flames with either pure liquid compounds or binary blends of these liquid compounds decreases monotonically with increasing carbon dioxide mole fraction X CO 2 in the oxidizer stream. Concomitantly, the flame length increases as the soot zone in the flame is shifted downstream. This effect can also be observed in Figs. 2

. At any toluene mole fraction, CO

2 dilution of the coflowing oxidizer results in a substantial decrease in integrated peak soot volume fraction in these diffusion flames. The other significant impact of CO 2 dilution is the downstream shift of the soot inception region, therefore of the integrated peak soot volume fraction, which results in longer flames.

As discussed in detail in Ref. 24, the addition of CO 2 affects soot formation through three main paths, i.e., dilution, thermal, and chemical effects.

In the region where soot formation dominates soot oxidation, these three effects accumulate to delay the soot inception process, which leads to the aforementioned shift. In the region where soot oxidation dominates soot formation, the third effect induces an opposite tendency to the first two as it hastens soot oxidation slightly. However, this third effect does not significantly appear.

The reaction CO

2 + H ! CO + OH [19] is primarily responsible for the chemical effects of CO

Combined effects

As already mentioned in section 3.1, Figure 3 shows the monotonic and potentially non-linear trends of YSI exhibited by the binary mixtures investigated. Therefore, a model scaling the absolute YSI as a function of X T ol. and X CO 2 is likely to be relevant. It will provide with a quantitative model predicting single effects and combined effects of both parameters on the YSIs of the gasoline surrogate fuels investigated. This would be a valuable benchmark for the validation of CFD simulations together with a simple engineering tool. To this end, the second-order polynomial in Eq. ( 7)

is expected to decently model the non-linearity that can be observed in Fig. 3:

Y SI = a T ol. ✓ 1 X T ol. X T ol. ◆ 2 + a CO2 ✓ X CO 2 X CO 2 ◆ 2 + a T ol./CO2 ✓ 1 X T ol. X T ol. ◆ ✓ X CO 2 X CO 2 ◆ + b T ol. ✓ 1 X T ol. X T ol. ◆ + b CO2 ✓ X CO 2 X CO 2 ◆ + c (7)
where the variables X T ol. and X CO 2 have been reduced by their respective ranges investigated, i.e.

X T ol. =1 and X CO 2 =0.06, to equitably weight their fluctuation spans. As pure toluene (X T ol. = 1) is the common fuel investigated within the range of conditions covered with both binary mixtures, the second-order polynomial is derived as a function of 1 X T ol. . The coefficient c, set by the measurements for X T ol. = 1 and X CO 2 =0, is then expected to be equal to the index of toluene Y SI T ol. for both n-heptane/toluene and isooctane/toluene mixtures.

The coefficients in Eq. ( 7) are determined using a trust-region-reflective algorithm [36] implemented into the Matlab environment to minimize the least-squares between the results of Eq. ( 7) and the set of YSIs measured.

The coefficients are reported in Tab. 3. In addition, a model such as the one expressed in Eq. ( 7) enables confidence intervals (CI) to be stated [24].

In Fig. 3, the red and green bands represent the 95% confidence intervals for blends of n-heptane/toluene and isooctane/toluene, respectively. Their lower and upper bounds are evaluated using Eq. ( 8):

CI (X T ol. , X CO 2 ) = Y SI (X T ol. , X CO 2 ) ± 1.96 ⇥ max (8)
where Y SI is the mean value of the YSIs measured at a given set of conditions (X T ol. , X CO 2 ) that are averaged over the three alternatives for YSI evaluation and both X vap,i . While, max is the maximum value among the standard deviations computed at every set of conditions (X vap,i , X T ol. , X CO 2 )

investigated. The interval ±1.96 ⇥ max represents 95% of the area under the standard distribution curve [START_REF] Kenney | Confidence limits for the binomial parameter[END_REF]. Indeed, the deviation of every value from this fit is considered a contribution to the uncertainty in the YSI measurement. Therefore, the confidence intervals reported in Fig. 3 represent the uncertainties of the YSIs as derived from Eq. (7).

Except the mixture of n-heptane/toluene at X T ol. =0.2 and X CO 2 =3%, all the measured YSIs fall inside the zone encompassed between the upper and lower bounds of the confidence intervals. Thus, the relevance of the fitting model expressed by Eq. ( 7) is assessed.

Figures 4 and5 exhibit the surfaces resulting from Eq. ( 7) with the coefficients obtained for the mixtures of n-heptane/toluene and isooctane/toluene, respectively. A set of YSI isolines are shown. Practical use of this kind of representation could be relevant as it can provide with the information for the selection of a range of conditions in terms of fuel mixture (X T ol. ) and CO 2 dilution of air (X CO 2 ) that would leave the amount of soot produced below a requested threshold.

For the blends of n-heptane/toluene, the major influence of the linear trends reveals as the isolines are quite straight in Fig. 4. This is evidenced by the relative magnitudes of the coefficients reported in Tab. 3: |a T ol. /b T ol. | n heptane =0.56 and |a CO2 /b CO2 | n heptane =0.10. Nonetheless, the non-linear dependency of YSI on X T ol. slightly grows as X CO 2 increases. This can also be observed in Fig. 4. The upper isolines shift towards the higher toluene mole fractions at a higher rate than the lower isolines when X CO 2 increases. Indeed, a T ol./CO2 and b CO2 have opposite sign. Therefore, when

X CO 2 increases, the cross term ⇣ 1 X T ol. X T ol. ⌘ ⇣ X CO 2 X CO 2
⌘ weakens the linear term to a significant extent as |a T ol./CO2 /b T ol. | n heptane =0.27.

For the blends of isooctane/toluene, the substantial influence of the nonlinear dependency on the blending parameter X T ol. exhibits as some of the isolines are really curved in Fig. 5. This is also evidenced by the relative mag-nitudes of the coefficients reported in Tab. 3: |a T ol. /b T ol. | isooctane =23.2 and |a CO2 /bCO2| isooctane =0.07. The latter ratio shows that for the blends of isooctane/toluene, the dependency on the CO 2 dilution of the air stream remains very linear. Like for for the blends of n-heptane/toluene, a T ol./CO2 and bCO2 have opposite sign. Therefore, the cross term

⇣ 1 X T ol. X T ol. ⌘ ⇣ X CO 2 X CO 2
⌘ further weakens the linear term when X CO 2 increases. However, this combined effect is much stronger in this case as the ratio of combined term and linear term is ten times the same ratio in case of n-heptane i.e. |a T ol./CO2 /b T ol. | isooctane = 2.7 while |a T ol./CO2 /a T ol. | isooctane =0.12.

Discussion

In the present study, no discernible synergistic effect on the sooting tendencies is revealed for blends of n-heptane/toluene while one exhibited for blends of isooctane/toluene. Interestingly, Kashif et al. [24] also observed no significant synergistic effect for blends of n-heptane/isooctane as the sooting tendencies of these blends showed strongly linear evolutions with isooctane mole fraction. This finding clearly highlights that the crucial -therefore controlling-role of aromatic hydrocarbons is very sensitive to the chemical composition of the fuel blend. It would also be prudent to remind here as already discussed in introduction that toluene makes up an important per- This mitigates ring formation rates. In contrast, the PAH precursor species experience later oxidizing conditions in the co-flow configuration which increases ring formation rates. This difference especially exhibits at lower toluene mole fractions, then weakens at higher toluene mole fractions, i.e. when the additional amount of PAH precursor species produced by the decomposition of isooctane is decreased.

Although not quantitative, this similarity further extends the consistency of the YSI methodology on the Santoro's burner.

Conclusions

The sooting propensities of binary mixtures of n-heptane/toluene and isooctane/toluene were evaluated in terms of Yield Sooting Indices (YSIs) in doped methane diffusion flames under the influence of CO 2 dilution of the air stream. YSIs were inferred from two-dimensional fields of soot volume fraction that a Laser Extinction Method probing the whole flame delivered.

While YSI increases in a quite linear way with toluene mole fraction in the blend of n-heptane/toluene, a synergistic blending effect revealed for the mixtures of isooctane/toluene. Indeed, the higher concentrations of propargyl radicals and vinylacetylene that are produced along the decomposition of isooctane hastens naphthalene formation and its growth to heavier PAHs.

This results in a higher rate of the soot production increase with toluene mole fraction for the blends of isooctane/toluene. Interestingly, these behaviors exhibited here over the co-flow burner match those reported for both kinds of blends in a counter-flow burner. This similarity further extends the relevance of the YSI methodology on the Santoro's burner.

Following a regular trend, the sooting propensity of every blend investigated decreases monotonically with increasing carbon dioxide mole fraction in the air stream. However, the combination of fuel blending and CO 2 dilution of air may result in synergistic effects. It especially contributes to emphasize the non-linear blending effect exhibited by the n-heptane/toluene blends at higher CO 2 dilution levels. A quadratic equation was set to decently reproduce these trends. Indeed, the prediction of the synergistic behavior can be crucial as most of the surrogate fuels used in numerical simulations and engine testing contain toluene mole fraction within the range where this effect is not negligible. Furthermore, this kind of approach could contribute to a careful methodology for the formulation of commercial gasoline.

Further efforts should now be dedicated to the measurement of temperature fields. These would provide significant insights into the understanding of the physical and chemical processes that lead to the aforementioned effects. Yield Sooting Indices of binary mixtures of n-heptane/toluene (blue dots) and isooctane/toluene (black dots), from pure nheptane or pure isooctane (X T ol. =0) to pure toluene (X T ol. =1), and for different CO 2 mole fractions in the air stream. Black solid lines are the interpolated YSI trends following Eq. (7).

List of Figures

The shaded bands represent the 95% confidence intervals (CI). 35 4 Surface modeling the evolution of YSI as a function of toluene mole fraction (X T ol. ) in the blend of n-heptane/toluene and CO 2 mole fraction (X CO 2 ) in the coflowing air. . . . . . . . . . 36 5 Surface modeling the evolution of YSI as a function of toluene mole fraction (X T ol. ) in the blend of isooctane/toluene and CO 2 mole fraction (X CO 2 ) in the coflowing air. . . . . . . . . . 37 34

Figure 1

 1 Figure 1 exhibits the soot volume fraction fields measured in methane flames doped with the vapors of different blends of n-heptane/toluene and isooctane/toluene (X vap,2 = 2.47 10 2 ) for some toluene mole fractions X T ol. and three different CO 2 contents of the oxidizer. Benzene's YSI value being

  2 addition. OH and O species attack soot precursors and delay soot nucleation. OH radical is known to attack mature soot particles resulting in lower soot volume fractions. The consumption of H species near the burner tip might slightly reduce the formation of the first benzene rings from toluene. This might also delay soot nucleation. These postulates are supported by the results shown in Figs. 1 and available results in the literature [19].

1 33 2

 33 Soot volume fraction fields in methane diffusion flames doped with different blends of: (a) n-heptane/toluene; (b) isooctane/toluene. X T ol. and X CO 2 represent the toluene mole fraction in the binary blends of liquid fuels and the carbon dioxide mole fraction in the air stream, respectively. The vapor mole fraction of the blend in the fuel stream is X vap,2 = 2.47 10 2 . . Evolutions of the integrated soot volume fraction F v (z) for different CO 2 mole fractions (X CO 2 ) in the air stream and for different mole fractions of toluene (X T ol. ) in the binary blends of: (a) n-heptane/toluene; (b) isooctane/toluene. The vapor mole fraction of the blend in the fuel stream is X vap,2 = 2.47 10 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3

Figure 1 :Figure 2 :

 12 Figure 1: Soot volume fraction fields in methane diffusion flames doped with different blends of: (a) n-heptane/toluene; (b) isooctane/toluene. X T ol. and X CO2 represent the toluene mole fraction in the binary blends of liquid fuels and the carbon dioxide mole fraction in the air stream, respectively. The vapor mole fraction of the blend in the fuel stream is X vap,2 = 2.47 10 2 .

Figure 3 :

 3 Figure 3: Yield Sooting Indices of binary mixtures of n-heptane/toluene (blue dots) and isooctane/toluene (black dots), from pure n-heptane or pure isooctane (X T ol. =0) to pure toluene (X T ol. =1), and for different CO 2 mole fractions in the air stream. Black solid lines are the interpolated YSI trends following Eq.(7). The shaded bands represent the 95% confidence intervals (CI).

Figure 4 :

 4 Figure 4: Surface modeling the evolution of YSI as a function of toluene mole fraction (X T ol. ) in the blend of n-heptane/toluene and CO 2 mole fraction (X CO2 ) in the coflowing air.

Table 3 :

 3 Coefficients of the polynomial fitting the experimental YSIs (see Eq.(7)).

		n-heptane/toluene	isooctane/toluene
	a T ol.	-10.81	-21.99
	a CO2	0.96	0.58
	a T ol./CO2	5.27	2.56
	b T ol.	-19.48	-0.95
	b CO2	-9.22	-8.75
	c	45.3	45.3