Bingqing Xie

Yuemin Zh
email: yuemin.zhu@creatis.insa-lyon.fr

Pei Niu

Ting Su

Feng Yang

Lihui Wang

Pierre-Antoine Rodesch

Loic Boussel

Philippe Douek

Philippe Duvauchelle

Super-Energy-Resolution Material Decomposition for Spectral Photon-Counting CT Using Pixel-Wise Learning

Keywords: X-ray CT, material decomposition, photon-counting detector, super energy resolution

Spectral photon-counting CT offers novel potentialities to achieve quantitative decomposition of material components, in comparison with traditional energy-integrating CT or dual-energy CT. Nonetheless, achieving accurate material decomposition, especially for low-concentration materials, is still extremely challenging for current sCT, due to restricted energy resolution stemming from the trade-off between the number of energy bins and undesired factors such as quantum noise. We propose to improve material decomposition by introducing the notion of super-energy-resolution in sCT. The super-energyresolution material decomposition consists in learning the relationship between simulation and physical phantoms in image domain. To this end, a coupled dictionary learning method is utilized to learn such relationship in a pixel-wise way. The results on both physical phantoms and in vivo data showed that for the same decomposition method using lasso regularization, the proposed super-energy-resolution method achieves much higher decomposition accuracy and detection ability in contrast to traditional image-domain decomposition method using L1-norm regularization.

I. INTRODUCTION

Spectral photon-counting X-ray computed tomography (sCT) provides additional X-ray spectrum information about material components in an organ or tissue, due to its ability of separating photons into different energy bins by using photon-counting detector (PCD). This ability makes possible efficient material decomposition of which the objective is to quantitatively measure different materials encoded by a pixel.

It is a huge challenge to obtain high accuracy of material decomposition in clinical applications, especially for

The associate editor coordinating the review of this manuscript and approving it for publication was Kathiravan Srinivasan . low-concentration materials [START_REF] Jacobsen | Dual-energy CT: Lower limits of iodine detection and quantification[END_REF], [START_REF] Xie | Image-domain material decomposition in spectral photoncounting CT for medical applications[END_REF]. This is linked to the principle itself of sCT. With respect to conventional CT or dual-energy CT (DECT) to some extent, sCT divides the assayed energy range into discrete energy bins. An increased number of energy bins correlates with an increased spectral resolution, although this might simultaneously contribute to a reduced signal-noise ratio (SNR), due to a smaller number of photons being collected in each energy bin. Many aspects can influence the accuracy of material decomposition, including exposure time, energy bin width, acceptable flux rate, etc. Although there is complex interdependency between the parameters, we will focus on energy bin that is a basic parameter for spectral CT in contrast to conventional CT. In other words, energy resolution of sCT is limited by a trade-off VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ between the total number of energy bins and the image quality in each bin, which restricts the accuracy of material decomposition and ability of detecting low-concentration materials.

Existing material decomposition methods can be divided into three categories: projection-domain [START_REF] Schlomka | Experimental feasibility of multi-energy photon-counting K -edge imaging in pre-clinical computed tomography[END_REF]- [START_REF] Zhang | Tensor-based dictionary learning for spectral CT reconstruction[END_REF], imagedomain [START_REF] Zhang | Tensor-based dictionary learning for spectral CT reconstruction[END_REF]- [START_REF] Wu | Improved material decomposition with a two-step regularization for spectral CT[END_REF] and one-step material decomposition [START_REF] Barber | An algorithm for constrained one-step inversion of spectral CT data[END_REF]- [START_REF] Mory | Comparison of five one-step reconstruction algorithms for spectral CT[END_REF], according to the sequence of material decomposition and image reconstruction, i.e. before, after or during. Most decomposition methods of all the three ways try to improve performance of material decomposition by exploiting the high correlation between spectral and spatial features embedded in the multi-energy images. However, almost no work can be found in the literature to enhance energy resolution of sCT by software. Actually, for DECT, there is an imaging method that aims to enhance the visibility of contrast agent by reconstructing the image at one chosen computed energy, called virtual monochromatic image (mono-E or VMS) [START_REF] Iyama | The usefulness of dual-layer spectral computed tomography for myelography: A case report and review of the literature[END_REF]. The performance of mono-E from DECT is determined by the material decomposition that is processed only on actual images, which is noisy and depends on the chosen computed energy.

In this paper, we propose to introduce the notion of superenergy-resolution (SER) into sCT in order to improve the accuracy of material decomposition in image domain. The idea is to learn the relationship between spectral features of simulated data and those of physical data, then applying it to other physical data. To do it, theoretical sCT images containing abundant spectral energy information, designated as SER images, are simulated according to physical phantoms. The mapping between actual sCT images and SER images are learned by utilizing a coupled dictionary learning (CDL) method in a pixel-wise way. SER images synthesized based on the learned mapping are then used for material decomposition in image domain.

II. MODEL AND METHOD

Common models of image-domain material decomposition are firstly shown in this section. Then, SER for sCT is introduced and modeled. At last, the proposed SER-based material decomposition is described.

A. MODELS OF IMAGE-DOMAIN MATERIAL DECOMPOSITION

Tomography reconstruction, which is a preliminary step for image-domain material decomposition methods, is separately processed at each individual bin. Reconstruction at the i-th energy bin tries to achieve the spatial image, i.e. the linear attenuation coefficients (LAC) µ( x, i) at pixel x:

µ( x, i) = argmin µ D(C(s i (u)), µ( x, i)) + R, ( 1
)
where s i (u) represents the recorded photons along the u-th ray in the i-th energy bin, C the mapping function between s i (u) and LAC, D the data-fidelity term, and R the regularization. The detailed relationship between µ and s, i.e. C in equation ( 1), is described by:

si (u) = R n 0 (E)d i (E)e -L(u) µ( x,E)dl dE, (2)
where si (u) denotes the mean value of s i (u), n 0 (E) the incident photon fluence from X-ray source, L(u) the u-th ray, µ( x, E) the LAC at energy E, and d i (E) the detector absorption efficiency or detector response function that can give the probability of a photon of energy E to be measured within bin i [START_REF] Schlomka | Experimental feasibility of multi-energy photon-counting K -edge imaging in pre-clinical computed tomography[END_REF]. For X-ray CT images, noise level is negatively related to the square root of the total number of detected photons [START_REF] Leng | Noise reduction in spectral CT: Reducing dose and breaking the trade-off between image noise and energy bin selection[END_REF].

In other words, sCT suffers from stronger quantum noise for narrower bin.

Material decomposition is then performed on the reconstructed multi-energy images:

µ( x, i) = M α=1 μmα (i)ρ α ( x), i = 1, . . . , B, (3)
where subscript α denotes the α-th basis material, μmα (i) the calculated i-th energy-bin effective mass attenuation coefficient (MAC), B the number of energy bins, M the number of basis materials, and ρ α ( x) the mass density at pixel x.

We rewrite image-domain material decomposition into its matrix form as:

Y = MX + N, (4)
where Y ∈ R B×N P designates the reconstructed multi-energy spatial images of µ with N P designating the number of pixels or voxels, M ∈ R B×M the decomposition matrix of µ m , X ∈ R M ×N P the decomposed basis material images of mass densities ρ and N the noise. Consequently the aforementioned trade-off between the number of energy bins and noise level leads to limited energy resolution in sCT. In practice, the number of energy bins of sCT is generally scheduled by manufacturers, e.g. 5 bins for Phillips sCT prototype [START_REF] Si-Mohamed | Multicolour imaging with spectral photon-counting CT: A phantom study[END_REF].

B. SUPER-ENERGY-RESOLUTION IMAGING IN sCT

We propose to improve the energy resolution by means of fully exploiting both simulation data and physical data. To do it, the notion of SER is introduced in sCT.

SER images are virtual multi-energy images with much more bins but noise-reduced compared to the actual sCT images. Fig. 1 illustrates the relationship between actual sCT, denoted by low-energy-resolution (LER) imaging, and SER imaging by respectively plotting the multi-energy values of corresponding pixels in the LER and SER images. For one pixel in the reconstructed LER images (e.g. yellow asterisks in the left of Fig. 1, designated as y L ∈ R B L ), it has 5 values for 5 energy bins, i.e. B L = 5. While for the corresponding SER pixel (yellow points in the right of Fig. 1, designated as y S ∈ R B S ), it can have much more values, e.g. 50 values for B S = 50 (30 -80 keV). Note that the choice of 50 bins corresponds to the most precise µ m we can get from NIST [START_REF] Saloman | X-ray attenuation coefficients (total cross sections): Comparison of the experimental data base with the recommended values of Henke and the theoretical values of scofield for energies between 0. 1-100 keV[END_REF]. That is, µ is sampled every 1 keV. The proposed SER imaging is achieved by introducing the assumption that if containing the same basis materials and concentrations, the mapping between LER and SER images is determined and can be learned. We realize it by learning the relationship between physical data and simulation data and then applying it to other physical data, as illustrated in Fig. 2. The physical data are measured and reconstructed as LER images, while the corresponding SER versions are produced by simulation with similar shapes and components but more bins. Then, LER and SER spectral features are extracted and their relationship will be learned by a mapping which can be utilized to synthesize SER images for other in vivo data.

C. SER IMAGE SYNTHESIS BASED ON COUPLED DICTIONARY LEARNING

We learn the mapping between LER and SER spectral features by utilizing CDL method in a pixel-wise way.

1) PIXEL-WISE TRAINING STRATEGY

It is inefficient to directly learn the mapping from the whole LER and SER images or patches. In general, learning morphological features of certain organ is another challenge that demands abundant images with various spatial information for training. It will be expensive to produce these kinds of physical data currently. Meanwhile, complex material components can make the learning more difficult.

We solve it by training the mapping from LER and SER pixel pairs, rather than the whole images or patches. With this strategy, the training samples only contain spectral information. Correspondingly, the physical phantoms used in our work can have simple shapes, but abundant material components.

Meanwhile, to cope with the difficulty of extracting features from noisy LER pixels, we enhance the learning by adding conventional CT images. The conventional CT images of physical phantoms are measured on the same sCT system but reconstructed with full spectrum (i.e. a single energy reconstruction), thus having lower noise level compared to that with single energy bin. The conventional CT and sCT pixels are then concatenated together as one LER sample for training. As a result, now y L ∈ R (B L +1) is a combined LER pixel.

2) LEARN THE MAPPING BY CDL

CDL was used in the field of super spatial resolution [START_REF] Yang | Image super-resolution via sparse representation[END_REF], [START_REF] Wang | Semi-coupled dictionary learning with applications to image super-resolution and photo-sketch synthesis[END_REF]. In the theory of dictionary learning, features of a pixel y can be represented as a sparse combination of atoms inside an over-complete dictionary D. That is, y = Dα, where α is the sparse codes with very few nonzero entries:

argmin α 0 , s.t. y -Dα 2 2 ≤ , (5)
where α 0 designates L 0 norm of α. Generally, solving optimization problem in ( 5) is NP-hard, and it can be instead solved using L 1 norm, as:

argmin α y -Dα 2 2 + λ α 1 , (6)
where λ denotes the Lagrange multiplier.

LER and SER pixel sets can then be separately represented using a LER dictionary D L and a SER dictionary D S :

argmin L Y L -D L L 2 F + λ L L 1 , argmin S Y S -D S S 2 F + λ S S 1 , (7)
where

Y L ∈ R (B L +1)×N t and Y S ∈ R B S ×N
argmin , D Y -D 2 F + λ 1 , s.t. Y = [ Y S Y L ], D = [ D S D L ], = [ S L ]. (8)
Meanwhile, a transform matrix is utilized to decrease the impacts of dramatic gap between simulation and physical data:

S = W L , (9)
where W denotes the transform matrix that was firstly introduced in super spatial resolution for stable cross-style image synthesis [START_REF] Wang | Semi-coupled dictionary learning with applications to image super-resolution and photo-sketch synthesis[END_REF]. The enhanced training model then becomes:

argmin D, ,W Y S -D S S 2 F + Y L -D L L 2 F + γ S -W L 2 F + λ W W 2 F + λ S S 1 + λ L L 1 . (10)
We solve the above model by alternating-directionmethod-of-multipliers (ADMM), which alternately updates It should be noted that atoms in the dictionaries learned using a pixel-wise strategy are energy-dependent, in contrary to spatial features in traditional super-spatial-resolution applications. Figure 3 illustrates the difference between super spatial resolution and super energy resolution by showing the learned features in practice. Atoms learned from a spatial image are morphological features, while for SER, atoms in the dictionary are energy-dependent and look more like µ curves.

After training, we can readily synthesize SER pixel ŷS from its LER version y L with the learned coupled dictionaries and transform matrix, by solving the following optimization:

SER(y L → ŷS ) : argmin ŷS ŷS -D S α S 2 F + y L -D L α L 2 F + γ α S -W α L 2 F + λ S α S 1 + λ L α L 1 . (11)
The above model is also solved by ADMM.

D. MATERIAL DECOMPOSITION USING SER IMAGES

The synthesized SER images are then utilized for material decomposition in image domain. For existing sCT, decomposition matrix M is generally initialized by the aforementioned μm [START_REF] Le | Least squares parameter estimation methods for material decomposition with energy discriminating detectors[END_REF] according to:

μmα (i) = E∈E i n 0 (E)d i (E)µ mα (E)dE E∈E i n 0 (E)d i (E)dE , i = 1, . . . , B, (12)
where µ mα (E) designates the theoretical MAC value and E i the width of the i-th bin. Note that μm (i) is an estimate of the true value of µ m (i) corresponding to the i-th bin. As a result, thinner bins can produce more precise μmα (i), i.e. more precise M. A more precise decomposition matrix suitable for the SER images is then needed. For LER, µ m of each basis material (one column of M in equation ( 4 decomposition using lasso as:

x = argmin x ŷS -M S x 2 F + λ x 1 . ( 13
)

III. EXPERIMENTS AND RESULTS

A. EXPERIMENTAL SETTING AND EVALUATION METRICS

The experiments were carried on real data from a Philips sCT prototype [START_REF] Van Der Werf | Coronary calcium scoring potential of large field-of-view spectral photon-counting CT: A phantom study[END_REF]- [START_REF] Si-Mohamed | Comparison of image quality between spectral photoncounting CT and dual-layer CT for the evaluation of lung nodules: A phantom study[END_REF]. This medical approved prototype has a 500mm FOV and a pixel pitch of 0.274mm at iso-center. Removable tubes containing different dilutions of iodine or gadolinium are inserted in the ROIs, as illustrated in Fig. 5b. The number refers to the concentration of each tube which covers from 0.1 mg/cc to 15 mg/cc. In our experiments, both iodine and gadolinium were measured for training and testing. Meanwhile, we use another in vivo dataset to further evaluate the proposed method. In our experiments, a rabbit was injected by gadolinium dilution, the image of which 5 and 15 minutes after the injection are respectively given in Fig. 6. For quantitative evaluation, tubes with known concentrations of gadolinium (2, 5, 10 mg/cc) are also placed beside the rabbit. The image of 5 min after injection is a little blurred due to movement.

All the spatial images in physical datasets were separately reconstructed at each energy bin via a common reconstruction method: conjugate gradient method (with Reconstruction Toolkit (RTK) [START_REF] Rit | The reconstruction toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the insight toolkit (ITK)[END_REF]). Note that the decomposition method of Philips prototype is projection-domain based, while our proposed SER-based method is performed in image domain. What's more, our reconstruction software (RTK) is also different from that used by Philips. In other words, we apply our own algorithms on prepped photon count data from the scanner, whereas other correction steps from the Philips image chain (e.g. proper handling of bad detectors) are not used. Hence the results will be different from the Philips' reconstruction results. The detector bin response function d i (E) in model ( 2) was provided by Philips. All phantoms were reconstructed into 9 slices and the size of each slice is 900*900. We learned the coupled dictionaries composed of 1024 (N D ) atoms from 400K (N t ) random LER-SER pixel pairs. On a machine equipped with Nvidia Tesla V100 and Intel I9 CPU, our method took 18h for training and 5 minutes for testing of one slice.

Both SER images synthesis and material decomposition were quantitatively evaluated. Normalized Root Mean Square Error (NRMSE) (i.e. NRMSE = RMSE/x for the image x) and peak signal-to-noise ratio (PSNR) were utilized to assess the quality of synthesized SER images. Note that NRMSE and PSNR were respectively calculated on the whole images of iodine and gadolinium, in order to assess the total synthesis performance of all concentrations. The smaller the NRMSE and the higher the PSNR, the more accurate the synthesis.

Mean and standard error of the mean (SEM) were calculated for each dilution to assess decomposition accuracy. Meanwhile, linear regression between the nominal (y-axial) and measured concentrations (using the calculated mean values, x-axial) was performed for all dilutions, as: y = ax + b, where a and b represent the slope and y-intercept, respectively. To further assess the detection ability of contrast agents, we also calculated the limit of detection (LOD) for both iodine and gadolinium. The LOD was described by the Clinical and Laboratory Standards Institute [START_REF] Moretti | CLSI EP17-A protocol: A useful tool for better understanding the low end performance of total prostate-specific antigen assays[END_REF] as:

LOB = xblank + 1.645δ blank , LOD = LOB -b a , (14)
where xblank denotes mean water signal detected in the images of decomposed contrast agents, δ blank the standard deviation of the measurement, and LOB the limit of blank. LOB denotes the highest apparent contrast agent concentration expected to be found when replicates of a sample contain no contrast agent, which was calculated using tubes of water in each phantom. LOD is taken as the lowest concentration of a contrast agent in a sample that can be detected.

B. TRAINING

NRMSE and PSNR (dB) of synthesized SER images at each energy bin are illustrated in Fig. 7. SER images of iodine have relatively higher accuracy than gadolinium, e.g. smaller NRMSE and higher PSNR, although both of them have good synthesis performance. Note that low-energy bins have relatively larger NRMSE and smaller PSNR, especially for bins of k-edge (i.e. energy bins containing k-edge). Also, we need to recover more details in the first two bins: 30 -50 keV and 51 -61 keV. The material decomposition results based on synthesized SER images are shown in Fig. 8. Visually, the proposed SER-based decomposition presents good detection ability according to the edge-preserving performance in training. It can substantially preserve edge information even for low-concentration iodine and gadolinium e.g. 0.5 mg/cc, as illustrated in the ROIs. Note that all the ROIs are shown with grayscale images in order to highlight morphological information.

To quantitatively assess the decomposition, linear regressions of iodine and gadolinium for all the concentrations are illustrated in Fig. 9. Note that there is no additional scatter correction for all the results of linear regression in this work. Both iodine and gadolinium have clearly high accuracy with respect to theoretical concentrations. The correlations between the measured and prepared concentrations are strongly linear for all dilutions (all R 2 ≥ 0.95). The slopes are close to 1 for both iodine (a = 0.96) and gadolinium (a = 0.97). The y-intercept values are very low, e.g. b = -0.26 for iodine and -0.38 for gadolinium. Additionally, the decomposition is very reliable because both iodine and gadolinium are always underestimated and their SEMs are small.

Table 1 lists LOD for both iodine and gadolinium, to further assess the detection limit. Both iodine and gadolinium can be detected for concentrations lower than 0.5 mg/cc in training. More specifically, the proposed SER method has better detection ability for iodine than for gadolinium.

C. TESTING

The LER-based material decomposition method using the same lasso regularization was compared to the proposed SER-based method. Both physical phantoms with more noise and artifacts and in vivo data were tested. SER images were through equation [START_REF] Wu | Improved material decomposition with a two-step regularization for spectral CT[END_REF]. Fig. 10 shows the NRMSE and PSNR (dB) of synthesized SER images at each energy bin for testing. As in the case of training, SER images of iodine have relatively higher accuracy than gadolinium, e.g. smaller NRMSE and higher PSNR. Compared to training, although testing has relatively higher NRMSE and smaller RSNR, the difference is small.

The material decomposition results of physical phantoms with more noise and artifacts are shown in Fig. 11. In contrast to LER method, the proposed SER method shows better performance in quantification and detection compared to LER method. Firstly, the proposed SER method gives much better morphological accuracy in terms of edge-preserving performance and much less noise, even for low-concentration materials, e.g. 0.5 mg/cc for both iodine and gadolinium. The edges are substantially preserved by SER method. While for LER method, the selected areas of low-concentration materials are overwhelmed by noise and cannot be detected.

Secondly, more quantitatively as shown in Fig. 12, linear regressions of both I and Gd illustrate that SER-based method has much smaller errors. SER method presents more precise slope (for iodine: a LER = 0.71 and a SER = 0.87; for gadolinium: a LER = 1.24 and a SER = 0.94) and y-intercept value (for iodine: b LER = 1.71 and b SER = -0.08; for gadolinium: b LER = 0.92 and b SER = -0.51). Meanwhile, SER method is more reliable with smaller SEMs and almost all decomposed values always exhibit the same underestimation trend. As can be observed in Table 1 that lists LOD values of different methods for both iodine and gadolinium, our SER method has much stronger detection ability for both iodine and gadolinium compared to LER method by producing smaller LOD values in testing. It has better detection ability for iodine than for gadolinium, as in the case of training.

Meanwhile, the decomposition results of different methods on in vivo data are given in Fig. 13. We identify calcium (Ca) inside spine as the basis material iodine in terms of their similar MAC under our energy-thresholds arrangement. Clearly, LER method has less accurate separation between gadolinium and iodine and has more noise for the two basis materials. In contrast, the SER method has better detection ability for both iodine (in the spine) and gadolinium (in tubes and organs). Mean and SEM of gadolinium in the tubes are listed in Table 2. Obviously, SER method leads to more accurate decomposition for gadolinium of concentrations covering from 10 mg/cc to 2 mg/cc.

IV. DISCUSSION

Our main goal is to improve the energy resolution of sCT in order to improve the k-edge quantification accuracy and detection ability. To achieve that we introduced the notion of SER in sCT to achieve SER imaging, by fully exploiting abundant spectral features embedded in both simulation and physical data. Preliminary results revealed that the proposed SER method enables accurate decomposition and high detection ability even for low-concentration materials.

The SER-based material decomposition method is reliable. Firstly, even for noisier conditions in training, the proposed SER method still enables high decomposition accuracy and detection ability. In our experiments, we trained the dictionaries on the phantoms with less noise and artifacts, but tested on the ones with stronger noise and more artifacts. However, the results still show better decomposition accuracy and detection ability of the proposed method compared to the LER method, especially for iodine. Secondly, accurate decomposition of in vivo data polluted by movement (5 min after the injection) also shows the reliability of our method. Moreover, the decomposition results of both iodine and gadolinium depict consistency, albeit being underestimated for almost all concentrations.

In this work, we did not consider the mixtures of multiple elements in the experiments. Nevertheless, materials detection mainly relies on the subsequent decomposition method if the input is SER images. In this sense, accurate SER images are primordial for high material decomposition accuracy. At the same time, the ability of the decomposition method with L1 norm to detect mixtures has been demonstrated in [START_REF] Xie | ROI-wise material decomposition in spectral photon-counting CT[END_REF], [START_REF] Xie | Material decomposition in X-ray spectral CT using multiple constraints in image domain[END_REF]. Hence, our method would work for the mixtures of interested elements (I or Gd). For other possible mixtures or multiple different contrast agents, e.g. bone/calcium basis, soft-tissue/water basis, corresponding training is necessary for ensuring the decomposition performance, which deserves further research in the future work.

Concerning the partial volume effect, which arises when more than one material type occurs in a voxel (or pixel) having a limited spatial resolution, it is a very interesting and pertinent issue in sCT. In our method, material decomposition is processed pixel by pixel, which implies that only material category and density are accounted for and spatial neighborhood is not involved. As a result, partial volume effect The interest of projection-domain and one-step methods lies in the fact that they generate directly material images from projection data without needing to first reconstruct spatial images. In contrast, image-domain decomposition has its own advantages of being able to both exploit spatial images and generate material images. On the other hand, projection data is not always available or not practical to use in clinical practice because of commercial confidentiality or practical constraints. In short, image-domain methods are hardware-independent and thus provide much flexibility than the other 2 methods.

Note that the accuracy of material decomposition is also conditioned by the limitations of PCDs, such as K-escape, pulse pileup, and charge sharing effects. Some PCDs can remove all the photons, which are under the lowest threshold before any charge sharing corrections causing incorrectly binned photons. Similarly, a shift in the energy of the measured photons (when over 40 keV) is possible, e.g., 60 keV photon is measured around 70 keV and is then binned in the wrong energy bin. To take into account the above-mentioned issues, a more sophisticated simulation of common problems seen in PCDs might be interesting, which would allow us to get insights into how the limitations of PCDs affect the accuracy of the SER-based material decomposition method.

Each single synthesized SER image is a virtual monochromatic image (mono-E). The mono-E image form DECT is generally reconstructed from a pair of decomposed basis material images (e.g. a soft-tissue-like basis material and a contrast agent). The performance of mono-E is determined by the material decomposition that is processed only on actual images (LER images for sCT in the present study). In contrast, SER images are directly synthesized using spectral features learned from both simulation and physical data. As a result, our method could contribute to more accurate mono-E imaging.

V. CONCLUSION

We have proposed a SER imaging method to improve image-domain material decomposition in sCT. The method is based on generating SER images from LER images. This is achieved through learning mapping between simulation and physical data in a pixel-wise way and applying it to other physical data. The results on both physical phantoms and in vivo data showed that the SER-based material decomposition method resulted in a clearly more accurate material decomposition and a higher detection ability as compared to the LER-based method. In the future work, the performance of other image-domain decomposition methods using SER images will be investigated. It would be also interesting to compare image-domain decomposition with projectiondomain methods such as that used by Philips sCT prototype.

FIGURE 1 .

 1 FIGURE 1. Linear attenuation coefficients: LER (left) and SER (right). Curves with different colors denote different pixels; corresponding LER and SER pixels have the same color.

2 F

 2 t respectively designate the LER and SER pixel sets with N t pixels in the training set, D L ∈ R (B L +1)×N D and D S ∈ R B S ×N D the LER and SER dictionaries containing N D atoms, L and S the corresponding sparse codes (i.e. sets of α in equation (5) and (6)), and Y -D the data fidelity term determined by Frobenius-Norm • F . Then, the connection between LER and SER features is established by jointly training D L and D S , formulated as:

FIGURE 2 .

 2 FIGURE 2. Illustration of SER based sCT imaging.

FIGURE 3 .

 3 FIGURE 3. Features Learned with dictionary learning method for super spatial resolution and super energy resolution. Note that the dictionary for sCT images (lower) is the SER dictionary. coupled dictionaries D, sparse codes and transform matrix W .It should be noted that atoms in the dictionaries learned using a pixel-wise strategy are energy-dependent, in contrary to spatial features in traditional super-spatial-resolution applications. Figure3illustrates the difference between super spatial resolution and super energy resolution by showing the learned features in practice. Atoms learned from a spatial image are morphological features, while for SER, atoms in the dictionary are energy-dependent and look more like µ curves.After training, we can readily synthesize SER pixel ŷS from its LER version y L with the learned coupled dictionaries and transform matrix, by solving the following optimization:

FIGURE 4 .

 4 FIGURE 4. Mass attenuation coefficients in M L and M S for iodine and gadolinium.

  ) designated as M L ) is an effective value at B L energy bins. In contrast, M for SER (designated as M S ) has much more energy bins, e.g. B S = 50. Both M L and M S have the same material components, e.g. iodine (I), gadolinium (Gd) and water in our experiment. Each column of M S is a more precise value directly retrieved from NIST, as illustrated in Fig.4.To assess the decomposition of the synthesized SER images, we choose a common image-domain decomposition algorithm using L 1 norm (lasso) regularization for both SER and LER images, which can reduce time complexity, preventing the model from over-fitting. Note that to illustrate the advantage of SER imaging, we focus on the comparison between SER and LER images, rather than different decomposition methods, i.e. only L1 method is considered in this work. Mathematically, we formulate the SER-based material 168488VOLUME 9, 2021

FIGURE 5 .

 5 FIGURE 5. (a) Reconstructed image (at 1st bin) of physical phantom for training (upper) and ROIs (areas inside red contour) for training and testing (lower). Note that the phantoms for training and testing have different noise levels due to the presence of extensional fat rings added to the testing phantom exterior. (b) inserts: labels designate the concentration in each disk (unit: mg/cc).

  The photon-counting detector has a 1848 × 64 grid. However we used a 2 × 8 binning mode to obtain 924×8 pixels with a size of 1 × 2mm. The experiments were performed at tube voltage 120 kVp and current 220 mA under rotation time of 1s with a small focal spot. The volumetric CT dose index (CDTIvol) was 22.2 mGy. The system acquires 2400 projections per rotation. It has 5 energy thresholds set to 30, 51, 62, 72, 81 keV. These values have been selected to optimize the Gadolinium k-edge imaging and are the same for all the physical data in the present work. This configuration leads to 5 energy bins: #1[30-51keV], #2[51-62keV], #3[62-72keV], #4[72-81keV] and #5[81-120keV].Reconstructed images from one single bin are shown in Fig.5a. Fig.5bdepicts the cylindrical inserts of the physical phantom. The upper image in Fig.5ais for training, while the lower two regions-of-interests (ROIs) of phantoms (i.e. all the areas inside red contours) are for training and testing, respectively. Training and testing sets do not share the same pixels. The reconstructed images for training and testing have different noise levels due to the presence of extensional fat rings added to the testing phantom exterior.

FIGURE 6 .

 6 FIGURE 6. The reconstructed image of rabbit after injection of Gd dilution (at 1st bin): (a) 5 minutes and (b) 15 minutes. Labels in ROIs designate the nominal concentration (mg/cc).

FIGURE 7 .

 7 FIGURE 7. The (a) NRMSE and (b) PSNR of synthesized SER images for I and Gd in training.

FIGURE 8 .

 8 FIGURE 8. The decomposition results of SER based method on physical phantoms in training. Note that ROIs are shown in grayscale to highlight morphological information. Labels in ROIs denote the dilution concentration (mg/cc).

FIGURE 9 .

 9 FIGURE 9. The linear-regression results of iodine and gadolinium for SER-based material decomposition in training: linear regression model, R-squared (R 2 ) and standard error of the mean (SEM).

FIGURE 10 .

 10 FIGURE 10. The NRMSE and PSNR of synthesized SER images for I and Gd in testing.

FIGURE 11 .

 11 FIGURE 11. (a) The decomposition results of different methods on physical phantoms in testing; (b) the ROIs in (a) for low-concentration materials.

FIGURE 12 .

 12 FIGURE 12. The linear regression results of (a) iodine and (b) gadolinium for different material decomposition methods in testing: linear regression model, R-squared (R 2 ) and standard error of the mean (SEM).

FIGURE 13 .

 13 FIGURE 13. The decomposition results of different methods on in vivo data in testing: (a) 5 min after the injection of contrast agents; (b) 15 min after. Note that calcium in spine is treated as iodine, and the basis material gadolinium exists in tubes and organs. The Numbers ('2', '5' and '10') refer to the concentrations in unit mg/cc.

TABLE 1 .

 1 The LOD of different material decomposition methods on physical phantoms.

TABLE 2 .

 2 The mean and standard error of the mean (SEM) of each tube in in vivo data (15 min after injection of gadolinium).

 VOLUME 9, 2021

VOLUME 9, 2021

 VOLUME 9, 2021

 VOLUME 9, 2021

ACKNOWLEDGMENT

The authors thank Dr. Klaus Erhard and Dr. Yoad Yagil for the careful reading of the manuscript and many helpful comments.

This work was supported in part by the Program Hubert Curien Partnerships (PHC)-Cai Yuanpei 2018 under Grant 41400TC, in part by the International Research Project Modélisation et traitement d'images et du signal pour la Santé (METISLAB), in part by the Committee of Science and Technology of Shanghai under Grant 19510711200, in part by the European Union Horizon 2020 under Grant 643694, and in part by the Project Point Accueil Installation (PAI) Region Auvergne-Rhone-Alpes under Grant 2000688501-40890.

PHILIPPE DUVAUCHELLE is currently an Associate Professor at LVA, INSA Lyon. His general skills concerns physics and a NDT methods. He is particularly an expert in techniques using X-ray, from simulation and modeling of X-rays-matter interactions up to images and signals analysis, through various types of conventional and nonconventional X-ray imaging. VOLUME 9, 2021


domain: phys sdu

text:
Emilie Hardy

Agnes Levy

Manuel Rodrigues

Pierre Touboul

Gilles Metris

G Métris

SOFTWARE TOOLS FOR THE VALIDATION OF THE IN-FLIGHT CALIBRATION PERFORMANCE OF THE MICROSCOPE SPACE MISSION

Keywords: MICROSCOPE, test of the Equivalence Principle, in-orbit calibration, measurement processing

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The MICROSCOPE space mission aims at testing the Equivalence Principle (EP). This principle is at the basis of General Relativity and has for main consequence the Universality of Free Fall: the acceleration of an object in a gravitational field is independent of its mass and its internal composition. It leads to the equivalence between the inertial mass (which measures the resistance of an object to a modification of its motion) and the gravitational mass (coupled to the gravitational field). The Universality of Free Fall has been tested throughout the centuries with an improving accuracy. Recently, sophisticated torsion-balances or Lunar laser ranging experiments have been conducted, leading to a record accuracy of a few 10 -13 (Gundlach et al. (2009)). However, the accuracy of these on-ground experiments is limited by the terrestrial environment's numerous perturbations. Some unification theories, representing alternatives to General Relativity and aiming at merging gravitation with the three other fundamental interactions, expect a violation of the EP at a level never tested before (Damour et al. (2002)). Being performed in space, the MICROSCOPE mission will be able to overcome the terrestrial perturbations in order to test the Equivalence Principle with an accuracy of 10 -15 (Touboul et al. ( 2001)).

MICROSCOPE is a microsatellite developed by CNES with a contribution of ESA. It will orbit around the Earth at low altitude for a two years' mission. The launch of the satellite is scheduled for 2016. The on-board payload is composed of two differential electrostatic accelerometers developed by ONERA, each one consisting in two concentric cylindrical test masses in electrostatic levitation. The position of the masses is detected using a capacitive measurement, and a control loop keeps them at the center of the accelerometer cage. Both masses are submitted to the same gravitational field. For one of the differential accelerometer, the two masses are made of different materials. A difference between the measured forces applied to maintain them on the same trajectory would therefore indicate a violation of the Universality of Free Fall, and thus of the Equivalence SF2A 2013

Principle. The second accelerometer includes two test masses of the same composition and enables us to test the measurement accuracy.

The accuracy of the measurement exploited for the EP test is limited by our a priori knowledge of the instrument physical parameters. It is necessary to finely characterize them in order to correct the measurement. An in-orbit calibration of the instrument is therefore needed. Appropriate calibration procedures have been determined and their analytical performances evaluated. These procedures have been implemented in a numerical simulator in order to validate them numerically (section 2). The estimation of the instrumental parameters is used to correct the scientific measurement and therefore to evaluate the EP violation parameter with the best accuracy. A mission simulator has been developed in order to simulate the measurement processing and thus to validate the mission performances (section 3).

2 In-orbit calibration of the instrument

Necessity and principle of the in-orbit calibration

The ideal measurement of the MICROSCOPE instrument corresponds to the electrostatic acceleration applied to the internal or the external test mass to keep it relatively motionless at the center of the accelerometer. The real measurement provided by the instrument differs from the ideal measurement because of instrumental parameters: scale factors, couplings, misalignments and quadratic parameters.

The measurements provided for the internal and the external test masses are combined in a common mode (the half sum of both measurements) and a differential mode (the half difference of both measurements). The common mode acceleration is used as a command for the satellite Attitude and Orbit Control System (AOCS). This system compensates for the external non-gravitationnal perturbations, such as the residual drag or the Earth and Sun radiation pressure. The differential mode is used for the EP test.

This differential acceleration includes the potential violation signal due to a non-null EP violation parameter

δ EP = mg2 mI2 - mg1 mI1
, where m g is the gravitational mass and m I is the inertial mass respectively of the internal test mass (index 1) and external test mass (index 2). A second order model of the differential measurement has been defined [START_REF] Hardy | Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics[END_REF]. The form of the acceleration measurement is:

Γ mes,dx = 1 2 • g x • δ EP + Γ inst,x (2.1)
with:

• g x the component of the gravitational acceleration along the X axis of the instrument;

• Γ inst the perturbative acceleration reducing the accuracy of the EP parameter determination. This acceleration is composed of a stochastic part (noise of the instrument and of the drag-free system) that can be reduced by integrating the measurement over a large duration, and an harmonic part that is due to instrumental parameters (scale factors, off-centering between the two test masses that introduces effect of the Earth gravity gradient and the inertia gradient, misalignements and couplings between the axes of the two sensors or between the axes of the differential accelerometer and the axes of the instrument, quadratical parameters that introduces non linear terms).

The instrumental parameters reduce the accuracy of the estimation of δ EP . The best performances possible for the instrument construction and integration in the satellite lead to a contribution of these terms incompatible with the mission's accuracy objective. It is therefore necessary to correct the measurement. An accurate estimation of the instrumental parameters is not possible on-ground, where the instrument is saturated. That is why in-orbit calibration sessions are needed. The in-orbit calibration of the MICROSCOPE instrument consists in amplifying the effect of the parameter to be estimated so that the corresponding term is predominant in the measurement, and thus easily measurable. For this purpose, a strong signal at a known frequency is used: either a naturally strong signal is used (the Earth gravity gradient), or a strong signal is created by oscillating the satellite (through the AOCS system) or the test masses (through their control loop) at a chosen frequency and along or around a chosen axis.

A specific calibration procedure has been determined for each instrumental parameter. Each calibration session will last 10 orbits in order to reduce the noise to an acceptable level, and the complete calibration will be run at least twice during the mission, so as to take into account a potential drift of the parameters value with time.

The performances of the calibration procedures have been evaluated analytically [START_REF] Hardy | Advances in Space Research Levy[END_REF]. The results are compatible with the specification. However, this analysis cannot take into account all the effects that impact the calibration result, in particular the dynamics of the systems. In order to validate their performances numerically, the calibration procedures have been implemented in a simulator.

The calibration simulator

The calibration simulator has been developed in Simulink (figure 1). It simulates the instrument, the satellite AOCS loop and the external environment. The simulated differential measurement is processed to estimate the parameter to be calibrated: a pass-band filter is applied around the calibration frequency and a least squares inversion is used. A statistical study based on 100 simulations with a random initialization of the parameters to be calibrated has been conducted. The resulting accuracy on the instrumental parameters estimations is compatible with the specifications [START_REF] Hardy | Advances in Space Research Levy[END_REF]. The performances of the foreseen calibration procedures are thus numerically validated. The in-orbit calibration provides an estimation of the instrumental parameters. These estimated values are SF2A 2013 used to correct the EP test measurement. In order to validate the complete mission scenario, including the inorbit calibration, the EP test session, the correction and the processing of the scientific measurement, a mission simulator has been developed in FORTRAN starting from a simulator of the OCA initially dedicated to the EP test sessions. The mission simulator enables us to simulate calibration sessions as well as EP test sessions, and the data processing for the estimation of the parameters (including the estimation of the instrumental parameters for the calibration, the correction of the scientific measurement using the estimated values of the instrumental parameters and the estimation of the EP violation parameter).

Measurement correction

In contrary to the calibration simulator, the mission simulator takes into account the orbital motion of the satellite. Therefore, an input of the simulator is an ephemerid file. This file provides the speed and acceleration of the satellite along its orbit. For each point of the orbit, the gravity acceleration and the Earth gravity gradient are computed from a model of the Earth gravity potential. Concurrently, the attitude of the satellite is computed. Three attitude laws are possible: inertial pointing of the satellite, spinning pointing, or oscillation of the satellite around one axis chosen for the calibration. These data are used to simulate the measurement of the electrostatic acceleration provided by the instrument. The simulated measurement is then processed by confronting it with a model of the measurement. The form of the model is:

Γ mes,dx = k a k • p k (3.1)
with:

• p k one of the parameters (instrumental parameters or EP violation parameter);

• a k its partial derivate.

The measurement is corrected using the best estimation pk of the parameter p k available for the analysis:

Γ mes,dx,corr = k a k • (p k -pk ) (3.2)
In practice, the estimated value of the instrumental parameters, known from a previous calibration, will be used. After correction, the measurement is processed with a least squares method. The form of the measurement model is:

Y = A • X (3.3)
with:

• Y the observation vector (the measurement);

• X the vector of the parameter(s) to be estimated;

• A the partial derivates matrix.

The least squares solution is:

X = (A T • W • A) -1 • A T • W • Y (3.4)
with W a symmetrical weight matrix. The variance of the estimator X is minimal if W is chosen so that W T •W is the inverse of the measurement noise's variance-covariance matrix. This solution is optimal, but induces some difficulties. The variance-covariance matrix is not known: a model of the noise is therefore necessary. Moreover, this matrix has a high dimension (equal to the number of observations), therefore the inversion is difficult. The solution chosen for this study consists in applying a Fourier transform on both sides of the equation before the inversion. In the spectral domain, only the frequencies around the calibration frequency or the EP test frequency are selected. Thus, the number of observations is significantly reduced. Moreover, the noise spectral density can be considered flat on these narrow frequency bands.

The mission simulator has been used to simulate the in-orbit calibration sessions followed by an EP test session, with a simulated EP violation parameter equal to 10 -15 . The EP test measurement has been corrected using the estimated values of the instrumental parameters, and the EP violation parameter has been estimated. The results of this estimation are gathered in table 1. Thanks to the calibration, the estimation accuracy after the measurement correction is better than 10 -15 , and thus compatible with the MICROSCOPE accuracy objective.

Fig. 1 .

 1 Fig. 1. The calibration simulator: the yellow blocks correspond to secondary inputs designed to oscillate the satellite or the test masses for the calibration.

Fig. 2 .

 2 Fig. 2. The mission simulator.

ONERA, FR-92322 Châtillon, France

Universit de Nice Sophia-Antipolis, Centre National de la Recherche Scientifique (UMR 6526), Observatoire de la Cte d'Azur, Goazur, 06130 Grasse, France c Société Francaise d'Astronomie et d'Astrophysique (SF2A) 2013

Conclusions

The budget of the measurement equation before calibration does not comply with the objective of the EP test accuracy. Several in-flight calibration sessions are therefore necessary during the space experiment in order to estimate the instrumental parameters that limit the measurement accuracy, so as to correct the measurement.

Parameters to be calibrated have been identified and appropriate methods of calibration have been proposed.

In order to validate the calibration process, a simulation software including models of the instrument and of the satellite drag-free system has been developed and the calibration methods have been implemented in the simulator. The results are compatible with the specifications.

For the data processing validation, a second software has been developed, enabling us to test the entire mission scenario, including the correction of the measurement with the parameters estimated during the calibration process. The simulated estimation of the EP parameter is compatible with the mission accuracy objective and validates the mission performances.

The simulators can also be used to study other effects that have not been taken into account in this paper. For example, missing data can occur because of teletransmission errors. The missing data can be simulated with the mission simulator in order to define its impact on the EP parameter estimation and to test appropriate measurement processing for the reduction of this impact.

This study prepares for the Scientific Mission Center (CMS) which is currently under development at ONERA and whose mission will be to validate the operational data and process the scientific data of the MICROSCOPE mission.

The authors wish to thank the MICROSCOPE teams at CNES and ZARM for the technical exchanges. This activity has received the financial support of ONERA and CNES.


Clara Rivero
email: crivero@mail.ucf.edu

Kathleen Richardson

Robert Stegeman

George Stegeman

Thierry Cardinal

Evelyne Fargin

Michel Couzi

Characterisation of the performance parameters of some new broadband glasses for Raman amplification

New glasses which exhibit a broadband Raman spectral response have been prepared and characterised for their Raman cross section and gain performance. Results on three different glass compositions in the system [(100-x)NaPO 3 -xNa 2 B 4 O 7 ]:TiO 2 /Nb 2 O 5 =1 are summarised. The impact on the Raman properties of changing the phosphorus to boron ratio while maintaining the TiO 2 /Nb 2 O 5 molar ratio equal to one, was examined. Polarised spontaneous Raman measurements are presented, and a comparison between the experimentally obtained Raman gain coefficient and that predicted from spontaneous Raman cross section measurements is made. Results obtained show that these glasses exhibit similar or slightly higher Raman gain coefficients than that of SiO 2 , with a spectral bandwidth almost five times greater (~1300 cm -1 ).

The demands for high speed data communication and the expected overload of data traffic requires the existing telecom window, currently operating in the C and L bands, to expand over a wider range of wavelengths. Such expansion to broader bandwidth would provide access to an additional number of channels allowing increased transmission of optical data along the optical fibre. This in turn requires an enhanced fibre transmission window, which until recently has been limited by the water absorption peak at around 1•4 µm. However, new improvements in the fibre optic industry have led to the advent of the AllWave fibre, the first full spectrum fibre designed for optical transmission over the entire telecom window of 1260-1625 nm. [START_REF] Refi | Optical fibres for optical networking[END_REF] Despite these advances, as information is transmitted along the optical fibre, the data signal still becomes depleted due to intrinsic losses in the fibre. Hence, an amplification medium is necessary to provide amplification of the signal in the network. Ideally, such a medium would amplify across the full spectral bandwidth of the signal, with minimal loss. Today, erbium doped fibre amplifiers (EDFA) are considered the primary means of amplification in existing optical networks; however, EDFAs are limited to amplification in the C and L telecom bands. [START_REF] Desurvire | Erbium-doped fibre amplifiers: principles and applications[END_REF] This leaves Raman amplifiers as the main medium available for amplification throughout the entire telecom window, given that the amplification bandwidth and position of the gain spectrum of a Raman amplifier is determined by the pump laser. [START_REF] Islam | Raman amplifiers for telecommunications[END_REF] Consequently, this enhancement to the telecom window has motivated the search for new materials as Raman amplification media. These materials should possess a broadband amplification window to allow minimisation of cost and components in the optical network. In other words, the broader the spectral response of these new materials, the more channels that can be amplified using a single amplifier and pump laser. Such a configuration would make the optical network more cost efficient and the Raman amplifier more attractive to future implementation in a real optical network. Of equal importance however, are crucial factors such as the trade off between the gain coefficient and the optical losses which must be taken into consideration to provide the best figure of merit available for a candidate material and its desired application.

Taking into consideration the different attributes defined above for new candidate gain media, we have examined a series of three different glass compositions in the system [(100-x)NaPO 3 -xNa 2 B 4 O 7 ]:TiO 2 /Nb 2 O 5 =1 in order to address and characterise their Raman gain response. These glasses have been selected because of their overall broadband spontaneous Raman spectral response. As will be shown, the substitution of boron for phosphorus and the ratio of these two components in the glass matrix allowed us to obtain a more uniform/flat spectral gain response that spans over 1300 cm -1 (~40 THz).

Experimental procedures

Glass making Glasses in the system [(100-x)NaPO 3 -xNa 2 B 4 O 7 ]:TiO 2 / Nb 2 O 5 =1, where x=5, 10 and 13 mol% (see Table 1) were prepared from high purity NaPO 3 (99•99% Aldrich), Na 2 B 4 O 7 (99•99% Aldrich), TiO 2 (99•995% Alfa Aesar), and Nb 2 O 5 (99•998% Cerac). The glasses were melted under an oxygen atmosphere, in platinum crucibles at a temperature of 1150°C, for 5 min. After melting, the glasses were quenched onto a preheated carbon plate, and annealed at a temperature 40°C below their glass transition temperature (T g ). Finally, the glasses were cut and optically polished.

Measurements

The thermal properties of the three different glass compositions were measured with a differential scanning calorimeter (DSC) at a rate of 10°C/min. The maximum temperature measured was 675°C, and the glass transition temperature was obtained from the inflection point of the DSC thermal curve (i.e. maximum of the derivative curve).

The density of the three different glass samples was measured by the Archimedes method in diethyl phthalate at room temperature (24°C). The accuracy of the measurement is within 0•02 g/cm [START_REF] Islam | Raman amplifiers for telecommunications[END_REF] .

The linear refractive index of the glasses was measured using the Brewster's angle method at a wavelength of 632•8 nm. The recorded reflection signal from a lockin amplifier was fitted using Fresnel equations. The experimental error was found to be ±0•02.

The linear absorption spectra of the different bulk samples were obtained using a Cary 500 Scan UV-Vis-NIR Spectrophotometer.

The polarised (VV and VH) spontaneous Raman spectra of the different bulk samples were measured using a micro-Raman setup, with an excitation wavelength of 514 nm, as described by Rivero et al. [START_REF] Rivero | Quantifying Raman gain coefficients in tellurite glasses[END_REF] The incoming polarised laser beam was focused onto the front polished surface of the sample via a 100× microscope objective, with a spatial resolution of about 2 µm. A polariser and quarter wave plate (l/4) combination were used to select the polarisation direction (vertical, V or horizontal, H) of the scattered light. A backscattered geometry was used to collect the Raman signal, which was spectrally analysed with a spectrometer and a CCD detector, with a typical resolution of about 6 cm -1 . The Rayleigh line was reduced with a holographic notch filter.

The Raman gain coefficient of the three different bulk samples, specified for the selected wavelengths (or frequency shifts), was measured using the experimental setup described by Stegeman et al. [START_REF] Stegeman | Tellurite glasses with peak absolute Raman gain coefficients up to 30 times that of fused silica[END_REF] A 10 Hz repetition rate Nd:YAG laser, with an excitation wavelength of 1064 nm, producing 50 ps pulses, was used as the pump wavelength for the gain measurement. The signal (probe) beam was provided via an OPG/OPA system, with a spectral output ranging from 780-2300 nm. Both, the pump and probe beams were spatially and temporally overlapped in the 1-2 mm thick bulk sample. After the sample, a monochromator was used to filter the pump from the probe wavelength. Following the monochromator, the 45° polarised signal beam (with respect to the pump beam) was split by a polarising beam splitter and directed towards two identical germanium detectors. The parallel component (VV) of the signal (with respect to the pump) was amplified given the exponential behaviour of the Raman gain process, and was used to measure 10% gain in amplitude, while the orthogonally polarised signal (VH) was used as a calibration reference. [START_REF] Stegeman | Tellurite glasses with peak absolute Raman gain coefficients up to 30 times that of fused silica[END_REF]

Results

The glass transition temperatures (T g ), crystallisation temperatures (T x ) (where applicable), densities (r), and refractive indices (n) of the three different compositions of the [(100-x)NaPO 3 -xNa 2 B 4 O 7 ]:TiO 2 /Nb 2 O 5 =1 family examined in this project are given in Table 1. In this study, we deliberately changed the phosphorus to boron ratio whilst keeping the TiO 2 /Nb 2 O 5 molar ratio equal to one.

Notice that, going from x=5-13% (i.e increasing the boron to phosphorus ratio), there is an increase in the glass transition temperature, and for x=10%, a crystallisation event is observed in the temperature range where the thermal analysis curves were measured (up to 675°C). No crystallisation event was observed for compositions x=5, 13% as measured over this temperature range.

Table 1 also shows that as the molar percent of boron to phosphorus ratio increases (i.e increasing x), there is an inflection of the densities and refractive indices of these glass samples. Substitution of boron to phosphorus from x=5-10% causes a decrease of the volumetric weight of the glass structure, and an increase in the refractive index; however, as we increase the boron molar percentage further to x=13%, there is a further increase of the density, and a small change in the refractive index.

In addition to these changes in the physical properties, Figure 1 illustrates the linear absorption coefficient of the three bulk samples. Firstly, one can see that these glasses are transparent throughout the entire telecom window, which makes them suitable for this type of application. Secondly, one can observe a blue shift in the absorption edge of the spectra, as shown in the inset graph, for increasing boron (x) content.

Figure 2(a) illustrates the spontaneous, polarised (VV and VH) Raman spectra of the three glass compositions. The VH Raman signal has been multiplied by a factor of three for comparative purposes. Since, in the majority of cases, the Raman spectrum of the glasses is highly polarisable, all the main Raman features can be obtained from the VV spectra. Due to this fact, our attention is mainly focused on the VV spectra. It is worth noticing that as the boron concentration increases (i.e x increases), the spontaneous Raman signal becomes more polarisation dependant, given by a decrease in the intensity of the VH Raman spectra. The spectral Raman gain coefficient was estimated using the equations described by Rivero et al [START_REF] Rivero | Quantifying Raman gain coefficients in tellurite glasses[END_REF] and Lines. [START_REF] Lines | Raman-gain estimates for high-gain optical fibres[END_REF] A fused silica (Suprasil) sample from Heraeus and a Schott SF6 glass were used as reference calibrations using the absolute values of Raman cross section for SiO 2 and SF6 reported by Pan et al. [START_REF] Pan | Raman scattering cross-section and nonlinear optical response of lead borate glasses[END_REF] The spontaneous Raman cross section values were obtained after correction for Fresnel losses and internal solid angle, as described by Kato & Takuma [START_REF] Kato | Absolute measurement of Raman scattering cross-sections of liquids[END_REF] and Galeener. [START_REF] Galeener | The relative Raman cross sections of vitreous SiO2, GeO2, B2O3, and P2O5[END_REF] Figure 2(b) illustrates the spectral Raman gain response of the three different glass samples. Note that the main contribution to the Raman gain curve is given by the Raman peak at 900 cm -1 . Also shown in Figure 2(b) is the experimental Raman gain coefficient obtained using the setup described by Stegeman et al. [START_REF] Stegeman | Tellurite glasses with peak absolute Raman gain coefficients up to 30 times that of fused silica[END_REF] From the figure one can see that both the estimated and experimentally obtained Raman gain curves show similar spectral features. All the measurements were performed relative to the values of the spontaneous Raman cross-section and gain coefficient of fused silica at the 440 cm -1 Raman vibration, given by Pan et al [START_REF] Pan | Raman scattering cross-section and nonlinear optical response of lead borate glasses[END_REF] and Stolen & Ippen, [START_REF] Stolen | Raman gain in glass optical waveguides[END_REF] respectively.

Discussion

From the results obtained in Table 1, we can clearly observe that as we go from x=10-13% (corresponding to a change in the boron to phosphorus ratio of 0•44:0•60), there is clear evidence of structural reorganisation taking place in the glass matrix. This reorganisation gives rise to a transformation in the thermal stability of the glass (resulting in the emergence of a measured T x for the x=10% sample), and an inflection point (minimum) in the density (±0•02 g/cm 3 ) for the same glass in the series. This reorganisation and corresponding property evolution results from a transition from a phosphate rich glass network to a borate rich network. Here these statements do not refer to any form or evidence of phase separation, but rather as an indication of which is the dominant network former in the glass. Similar results associated with such a transition have been observed in the glass system (100-x)NaPO 3 -xNa 2 B 4 O 7 , [START_REF] Ducel | Physical and chemical characterisations of sodium borophosphate glasses[END_REF] where changes in structural configurations with variation of the B/P network were verified via NMR. [START_REF] Ducel | 31 P MAS and 11 B NMR study of sodium rich borophosphate glasses[END_REF] The reorganisation in the glass structure has been attributed to the transformation of four coordinated boron to three coordinated boron. [START_REF] Ducel | 31 P MAS and 11 B NMR study of sodium rich borophosphate glasses[END_REF] In the present study, the three coordinated boron content abruptly increases in the region where the glass network goes from a phosphate rich to borate rich (i.e x=10-13%) network, where the boron is less likely to participate in the glass matrix as linkages to the P-O-P bonds. The change in the role of the boron in the network causes the inhomogeneous introduction of titanium and niobium into the glass network, as when present in low quantities, these heavy metal species function as modifiers. Here, the introduction of titanium and niobium into the glass is more compatible for the case of a phosphate dominant network, where the less constrained network can accommodate the large ionic size of other species. These structural changes will be described in the interpretation of the Raman spectra below. Also note that the x=10% composition is less stable (important for fiberisation) since the difference between its T x and T g is smaller, as compared the other two compositions.

Another essential property for this type of application is the absorption of these materials in the wavelength range of interest. Figure 1 illustrates the linear absorption curve of these three glass compositions. Note that these glasses are transparent throughout the entire telecom window. Also the inset figure shows that, as we increase the boron to phosphorus ratio (i.e increasing x), there is a blue shift in the absorption spectra. It is important to also mention that these three particular glass compositions in the family [(100-x)NaPO 3 -xNa 2 B 4 O 7 ]:TiO 2 /Nb 2 O 5 =1 were chosen due to the fact that their respective parent glass network constituent exhibited the lowest OH related absorption as described by Ducal et al. [START_REF] Ducel | Physical and chemical characterisations of sodium borophosphate glasses[END_REF] Such absorption would be detrimental to any amplifier applications.

Shown in Figure 2(a) is the polarised spontaneous Raman spectra of these three glass compositions. It can be seen that a detailed analysis of the spontaneous Raman spectra of the three glass compositions studied is very complex due to the multiple contributions of the multiple Raman active species in the glass network. The designation of the main spectral features in the spontaneous Raman spectra of borophosphate glasses in the system (100-x)NaPO 3 -xNa 2 B 4 O 7 have been fully described elsewhere. [START_REF] Ducel | Structural study of borophosphate glasses by Raman and infrared spectroscopy[END_REF] Additionally, the Raman spectra of two similar glass families NaPO 3 -Na 2 B 4 O 7 -TiO 2 and NaPO 3 -Na 2 B 4 O 7 -Nb 2 O 5 have been investigated by Cardinal et al. [START_REF] Cardinal | Raman scattering and XAFS study of optically nonlinear glasses of the TiO 2 -NaPO 3 -Na 2 B 4 O 7 system[END_REF][START_REF] Cardinal | Correlation between structural properties of Nb2O5-NaPO 3 -Na 2 B 4 O 7 glasses and nonlinear optical activities[END_REF] Using these data, we can identify the main Raman bands given in Figure 2(a). Specifically, the main phosphate vibrations observed at around 685, 1014, 1170, and 1250 cm -1 , have been attributed to the symmetric and asymmetric P-O-P chain vibrations, and symmetric and asymmetric stretching modes of PO 2 entities, respectively. Note that no strong boron Raman bands are noticeable in the spectra since their respective Raman vibrations are convolved within the strong Raman bands for the other glass constituents. In addition to these bands, one can clearly observe a strong Raman vibration at 900 cm -1 , which is attributed to the vibration of the nonbonding Nb-O short bond, found in iso-lated NbO 6 octahedral units. Convolved within this band is also the Ti-O short bond Raman band at around 930 cm -1 . Other notable Raman signatures can be seen for features at 640-650, 740, and 830 cm -1 , which correspond to the symmetric NbO 6 and/or TiO 6 octahedral vibrations, and distorted Ti-O-Ti, and Nb-O-Nb chain vibrations, respectively. The low frequency (<400 cm -1 ) vibrations are associated with the vibrations of the heavy metal cations (i.e Ti +4 and Nb +5 ). [START_REF] Miller | The intensity of Raman scattering in glasses contain-ing heavy metal oxides[END_REF] Figure 2(b) shows the Raman spectral gain curve calculation for the three bulk samples (black), determined from the absolute spontaneous, polarised Raman intensity measurements illustrated in Figure 2(a). Note however, that the estimation of the spontaneous Raman cross section is more complex and hence less accurate (±15%) in this glass system, as compared to the calculations presented in our previous work on tellurite glass systems. [START_REF] Rivero | Quantifying Raman gain coefficients in tellurite glasses[END_REF] We attribute this slight increase in error to the complexity of the glass system, where there are multiple contributions from the different glass constituents to a particular Raman band, all of the same order of magnitude, providing a continuous Raman bandwidth. This is in contract to earlier results [START_REF] Rivero | Quantifying Raman gain coefficients in tellurite glasses[END_REF] on tellurites where only one glass constituent (TeO 2 ) dominated the gain response and cross section could be predicted and measured. Consequently, the individual deconvolution of the Raman spectral features of each constituent is required to provide better agreement between the calculated spectral Raman gain response, from the spontaneous measurements, and the experimentally obtained Raman gain coefficients, due in large part to the incoherent nature of the spontaneous Raman process as compared to stimulated Raman. Such detailed deconvolution has not been performed here.

Nonetheless, from the estimated Raman gain curve response one can clearly observe that the main spontaneous Raman spectral features are still present in the Raman gain spectra, provided the appropriate Bose-Einstein correction factor for the spontaneous spectra is used. [START_REF] Stolen | Issues in Raman gain measurements[END_REF] It is also evident in Figure 2(b) that the greatest contribution to the gain coefficient is given by the Raman peak at 900 cm -1 , which is attributed to the isolated NbO 6 octahedral vibration. This was expected since previous work [START_REF] Cardinal | Raman scattering and XAFS study of optically nonlinear glasses of the TiO 2 -NaPO 3 -Na 2 B 4 O 7 system[END_REF][START_REF] Cardinal | Correlation between structural properties of Nb2O5-NaPO 3 -Na 2 B 4 O 7 glasses and nonlinear optical activities[END_REF][START_REF] Miller | The intensity of Raman scattering in glasses contain-ing heavy metal oxides[END_REF] has shown that d 0 metal ions such as Ti +4 , Nb +5 , and W +6 are highly polarisable. A small enhancement of the gain coefficient at this particular vibration, was observed for the glass composition with x=5% (i.e largest phosphorus to boron ratio), since, as previously discussed, a phosphate network is more acceptable to the introduction of the d 0 ions, giving rise to a slight increase in the number density of titanium and niobium species that can be accommodated within the glass network. One can also notice however, that the overall Raman gain response for this composition (x=5%) is less uniform over the spectral range shown, as characterised by the gain coefficient in the low (up to 1170 nm), peak (1180 nm), and high (>1180 nm) regimes. As previously noted, as the boron concentration in the glass is increased, the glass network structure depolymerises (due to the above mentioned boron structural re-arrangement), giving rise to a more uniform, flat Raman gain spectral re-sponse, which is more desirable for broadband amplification applications.

Finally, the experimentally obtained Raman gain coefficients are also illustrated in Figure 2(b), and they show good agreement with the spectral gain data calculated from the spontaneous cross section measurements; the same spectral features are seen (within the experimental errors of both techniques). A measurable discrepancy between the experimental and predicted gain spectrum is observed in the low wavelength region, most notably for the x=13% sample. No explanation for this observation can be presented at this point; however, there should not be any signs of resonance enhancement at the wavelengths where these measurements were performed, given that both calculations were performed far away from the absorption edge of the sample (Figure 1). Further analysis is being undertaken to account for this discrepancy between the two results. Nonetheless, the Raman gain curves shown in Figure 2(b) reveal that the gain response of these glasses is of the same order of magnitude as fused silica, with an improved spectral bandwidth of about 40 THz, an increase in bandwidth of more than five times as compared to SiO 2 .

Conclusions

This paper presents Raman gain performance results on glasses within the system [(100-x)NaPO 3 -xNa 2 B 4 O 7 ]:TiO 2 /Nb 2 O 5 =1, with x=5, 10 and 13%.

Results show that there is an inflection of the density and linear refractive index, as the glass composition moves from a phosphorus rich to a boron rich network. This network transformation also impacts the thermal stability of the glass, with the x=10% composition exhibiting the least stability as regards fiberisation. Analysis shows that a structural transition similar to that seen in glasses without Ti/Nb occurs which exhibits a slight decrease in glass stability for x=10% composition.

From the polarised spontaneous Raman data we were able to obtain qualitative information regarding the glass structure, and estimations of the spectral gain curve for these three different glass compositions, based on spontaneous cross section measurements. Our estimations show that these glasses exhibit an overall spectral gain of the same order of magnitude as fused silica, with an improved spectral bandwidth that spans, continuously, over 1000 cm -1 , corresponding to 50 THz of spectral bandwidth available for amplification applications (when pumped at 1064 nm), an improvement of almost five times as compared with SiO 2 .

Figure 1 Figure 2 (Figure 2 (

 122 Figure 1. Absorption spectra

Table 1 .

 1

	[(100-x)NaPO 3 -xNa 2 B 4 O 7 ]: Thermal analysis	r (g/cm 3 )	n±0•02
	TiO2/Nb2O5=1	Tg (°C) Tx (°C) ±0•02 g/cm 3 632•8 nm
	x=5%	433	-	2•72	1•54
	x=10%	459	590	2•68	1•56
	x=13%	468	-	2•75	1•55

Acknowledgment

This work was carried out with the support of a number of research, equipment, and educational grants, including NSF grants ECS-0123484, ECS-0225930, INT-0129235, NSF IGERT grant DGE-0114418, and NSF-CNRS grant 13050. The authors acknowledge the assistance and support of our collaborators at ICMCB / LPCM, and the School of Optics / CREOL, as well as the financial support offered by the UCF Graduate Studies office and the Student Government Association (SGA). Also, special thanks to Schott North America-Regional R & D for providing us with the SF6 glass sample.


Clémence Houzé

Albane Saintenoy

Léa Pannecoucke

Mathieu Le Coz

Charlotte Cazala

Chantal De Fouquet

Estimations of hydrodynamical parameters of Fontainebleau sands using GPR monitoring of water infiltrations

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Hydrodynamic properties of soils is a key knowledge for developing hydrogeological models to many applications such as following the evolution of pollution plumes (Pannecoucke et al., 2019). These properties can be defined using the water retention curve Ө(h) as well as the hydraulic conductivity curve K(h), representing respectively the relationship between the volume water content of the soil, Ө, or the hydraulic conductivity of the soil, K, as a function of suction h. These two curves are usually modeled using the [START_REF] Mualem | Hysteretical models for prediction of the hydraulic conductivity of unsaturated porous media[END_REF]) -van Genuchten (1980) formalism, valid on a wide variety of soils, using 5 parameters. These parameters can be determined by different laboratory techniques among which drainage measurements with water hanging columns (WHC), as described in [START_REF] Dane | Laboratory determination of water retention[END_REF]. Many field methods exist based on monitoring water infiltration in different configurations (Porchet et Lafferere, 1935;[START_REF] Coquet | Caractérisation hydrodynamique des sols à l'aide de l'infiltromètre à disques[END_REF]. Another way is to use pedo-transfer functions linking hydrodynamic parameters to physical characteristics of the soil such as the particle size distribution used by Rosetta (Zhang and Schaap, 2017).

Many geophysical methods are applied with success for insitu water flow determination (Robinson et al., 2008). Among these methods, the GPR is very sensitive to the aqueous phase of the medium and is commonly used to measure the water content of soils [START_REF] Huisman | Measuring soil water content with ground penetrating radar[END_REF]. One method, proposed by Saintenoy et al. (2008) and taken up by [START_REF] Léger | Hydrodynamic parameters of a sandy soil determined by ground penetrating radar inside a single ring infiltrometer[END_REF][START_REF] Léger | Hydrodynamic parameters of a sandy soil determined by ground-penetrating radar monitoring of Porchet infiltrations[END_REF] is based on monitoring an infiltration bulb during a single ring or a Porchet type water infiltration experiment in a shallow borehole. In this second case, the arrival times of the reflection of the waves at the surface of the bulb are inversed via an optimization tool using numerical models coupling the hydrodynamic code SWMS2D (Šimůnek et al., 1994) and the electromagnetic code gprMax (Warren and Giannopoulos, 2016).

Here, we test this hydrogeophysical method, which we call Porchet-GPR, comparing it to the other conventionally used experiments in order to demonstrate its capacity or not to replace them. This method has been applied to different sites of Fontainebleau sands. The Porchet-GPR Mualem -van Genuchten (MvG) parameters are compared with parameters obtained with i) the Rosetta3 pedotransfer function from particle size measurements and ii) WHC measurements. The MvG parameters are the residual water content Өr, the water content at saturation Өs, the hydraulic conductivity at saturation Ks, two parameters for adjusting the water retention curves, n and α, and an initial water content Өi. We are looking for the MvG parameter set which best explains the two types of data, i.e. VOL and TWT.

The Porchet-GPR method

The Porchet-GPR set-up can vary in the depth of the borehole, its radius, the water layer thickness maintained at the bottom of the borehole and the distance between the center of the borehole and the antennas. For a given geometry, we calculate using a Python3 application, around 20,000 models by varying the sets of MvG parameters in wide ranges (Өi=0.09 or 0.1, Өr=0.03, 0.36 ≤ Өs ≤ 0.44, 2 ≤ n ≤ 10, 1 ≤ α ≤ 10 1/m, 0.06 m, 0.06 ≤ Ks ≤ 0.65 cm/m, 0.06 min). For each set of parameters, the water content maps corresponding to each infiltration time step are calculated using the unsaturated flow numerical modeling code SWMS-2D. These water content maps are converted into relative dielectric permittivity maps using in the case of Fontainebleau sands, the CRIM formula [START_REF] Birchak | High dielectric constant microwave probes for sensing soil moisture[END_REF]. The GPR trace corresponding to each water content map is calculated by gprMax assuming a Ricker signal of frequency 1 GHz, and the given acquisition geometry (homogeneous medium in which a borehole of radius and depth is defined with a layer of water at its bottom). For each set of MvG parameters, the traces calculated at each infiltration time are concatenated to obtain a simulated radargram. On this radargram, the arrival times of the reflection on the infiltration bulb (TWT) are identified on each trace.

Field experiments

Experiments were deployed over three 100-square-meter stretches of Fontainebleau sandy soils, Auffargis, Bilboquet and Poligny. At each site, i) we sampled the soil with 250 cm 3 calibrated cylinders, ii) we made a series of auger holes, Riverside type, in which soil samples of any volume were taken, iii) we conducted Porchet-GPR experiments in each of the auger holes.

Once in the laboratory, all of the samples were immediately weighed before going to the 105 ° oven for 48 hours. They were weighed after the oven to determine their initial water content and their porosity for the 250 cm 3 samples. The soil cylinders are then dipped in water for an hour before being used for hanging water column measurements in drainage conditions. Finally, all of the calibrated or non-calibrated samples are analyzed with laser diffraction spectroscopy to determine their particle size distribution.

Results

Fig. 2 shows the particle size distribution for all investigated sites. A slightly coarser particle size is observed for the sands of the Poligny site compared to the particle size of the sands of the two other sites. The granulometric data make it possible to estimate the MvG parameters thanks to the Rosetta3 pedo-transfer function. The minimal variations in particle size have a notable effect on the parameters n, α and especially Ks which can vary from simple to double (Fig. 3). The pedo-transfer function is therefore sensitive to even small variations in the sand and silt content, while the clay content remains constant. Thelarger particle size at Poligny results in a slight increase of MvG parameter values estimated by Rosetta.

The hanging water suspended column experiments give an estimate of the water retention curves of the samples taken at the three study sites (Fig. 4). All the experiments were carried out in drainage. As with the particle size results, the curve obtained with the samples from the Bilboquet site is slightly above that obtained in Auffargis, itself above that obtained in Poligny, the place where the samples had the coarsest particle size. However, the general appearance is similar for the three curves, very flat at the point of inflection, which indicates an extremely rapid desaturation of the sample after the air-entry value is obtained. Adjusting the measurement points by the equation of van Genuchten (1980) gives parameters in Table 1.

Concerning the Porchet-GPR experiments, we compare the measurements on the Bilboquet and Poligny sites, on which the same acquisition geometry was used: a borehole 30 cm deep, radius 4 cm, and 5 cm of water kept at the bottom of the hole for at least 6 minutes. GPR trace is acquired every 5 s. The cumulative volumes of infiltrated water over time are in Fig. 5 and the TWT of the radar reflection on the infiltration bulb are in Fig. 6.

There are two types of data: two-way travel times (TWT) and volumes (VOL). We calculated the square root of squared errors (RMSE) between each data set and the results of 20000 models obtained by varying the hydrodynamic parameters. Fig. 7 is showing fitting results, depending on the data taken for the fitting, for one Porchet-GPR experiment carried out at the Bilboquet site. From different data sets, it is clear that modeling Porchet-GPR with Rosetta parameters do not fit field TWT nor VOL. Keeping all models giving a TWT fit less than 0.3 ns, results in a wide range in n value. Table 2 summarizes the parameters fitting field data acquired at the Bilboquet site when n varies in range [3,5].

Conclusions

Using the MvG parameters obtained from the Rosetta pedotransfer function we do not mangae to model the Porchet-GPR experimental results. From fitting both TWT and VOL data, we see that the n parameter is not constrained at all. Using a n value comparable to the Rosetta one, we obtain similar values for α but a Ks value 3 times smaller than the one from Rosetta. The α value with Rosetta and Porchet-GPR methods is double than the WHC one. It might be explained from hysterese effect in draining vs wetting conditions. Similarly, difference between hydrostatic (for WHC) and dynamic (for Porchet-GPR) conditions might explain difference in n estimations.

Adding GPR data to classical hydrological techniques illustrates the difficulties in hydrodynamical parameter estimations.

Figure 1 :

 1 Figure 1: The Porchet-GPR set-up The Porchet-GPR experimental set-up is described in Léger et al. (2015) and Fig. 1. The water infiltrated volumes (VOL) are measured every 30 s during the infiltration experiment that last about 6 minutes in sands. The two-way travel times (TWT) of the reflected wave going around the bulb are picked on the radargram acquired with static antennas placed on the surface.

Figure 2 :

 2 Figure 2: Synthesis of laser granulometric measurements done on field samples from the three sites

Figure 6 :

 6 Figure 6: TWT of GPR reflections on the infiltration bulb for Porchet-GPR experiments at Bilboquet and Poligny sites.

Figure 5 :

 5 Figure 5: Cumulative infiltrated volume of water during Porchet-GPR experiments at Bilboquet and Poligny sites

Table 1 :

 1 van Genuchten parameters from fitting water hanging column measurements

Acknowledgments

This study is part of Kri-Terres project, supported by the French National Radioactive Waste Management Agency (ANDRA) under the French "Investments for the Future" Program.

 Pannecoucke, L., Le Coz, M., Houzé, C., Saintenoy, A., Cazala, C., & de Fouquet, C. (2019). Impact of spatial variability in hydraulic parameters on plume migration within unsaturated surficial formations. Journal of hydrology, 574, 160-168. Porchet, M., & Laferrere, H. (1935)