Laurent JEANNEAU 1 , Richard D. ROWLAND 2 and Shreeram INAMDAR 2 1 . Géosciences Rennes, CNRS-Université de Rennes 1, France; 2 . Plant and Soil Sciences, University of Delaware, Newark, DE, USA

Discussions and message to carry back

• In the 12Ha catchment, the POM mainly derived from litter, wetland A and forest O horizons. Their relative contributions changed along the particles size continuum and from one event to another.

• At the scale of the headwater catchment, the sources of POM changed with an increasing proportion of stream bed-derived POM downstream. Going downstream seems to decrease the differences along the particles size continuum.

• The composition of POC changed along the particle size continuum: different fractions from different sources.

• Investigation at the molecular scale using THM-GC-MS seems to provide fruitful data to determine the sources of POM transferred along nested catchments.

Investigation at the molecular scale

Thermally assisted hydrolysis and methylation

Macromolecules lysis

Step 1: Thermochemolysis

Step 2: Gas chromatography

Separation

Step 3: Mass spectrometry 1. Where does this POM come from?

2. Is it homogeneous along the particle size continuum?

3. How its source may change down a continuum of catchment drainage locations?

4. Is there a seasonal effect?

  site • Piedmont region of Maryland, USA • Cold-temperate climate • Annual rainfall: 1007 (2014), 1165 (2015) mm • Elevations from 48 to 116 m (89 ± 14 m; mean ± SD) • Slopes from 0 to 35% (7.4 ± 5.0 %; mean ± SD) Nested catchment 12 and 79 Ha Sampling of deposited sediment during 4 storm events Sampling of potential end-members Outlet of the 12 Ha catchment PCA treatment based on the distribution of the type of molecules Individuals: End-members Variables: % of fatty acids (FA), carbohydrates (CAR), phenols (PHE) among the analyzed compounds and the sum of the area per gram of analyzed sample. 2. Second PCA treatment based on the distribution of plant-derived markers Individuals: End-members Variables: Relative proportion of the eight most abundant PHE and plantderived FA. potential end-members. Specific areas are increased by spatial heterogeneity. The composition of POM changed along: (1) the catchment (differentiation between the samples from the 12 Ha and the 79 Ha catchments), (2) the size continuum from coarse to fine fractions and (3) the year. The 12 Ha catchment POM derived mainly from a combination of Litter, Forest O and Wetland A horizons. The contribution of stream bed increased downstream. (POM) plays a key-role toward biodiversity, ecotoxicology and within the C cycle. POM exported within the fluvial suspended matter load during storm events may represent up to 80% of its annual load.

Acknowledgments

This work benefits from NSF and USDA funding (grant numbers IIA-1330238 and 2015-37020-23585, respectively).


Marc Allain
email: marc.allain@fresnel.fr.

Jérôme Idier
email: jerome.idier@irccyn.ec-nantes.fr.

Efficient binary reconstruction for non destructive evaluation using gammagraphy

Keywords: numbers: 81.70.-q Inverse Problems Non destructive evaluation, tomography, binary reconstruction

come

Efficient binary reconstruction for non destructive evaluation using gammagraphy

Introduction

Reconstruction of binary volumes from noisy projections is an attractive tool for non destructive evaluation (NDE) of metal structures [START_REF] Chalmond | Tomographic reconstruction from non-calibrated noisy projections in non-destructive evaluation[END_REF][START_REF] Retraint | Three-dimensional regularized binary image reconstruction from three two-dimensional projections using a randomized ICM algorithm[END_REF]. It amounts to the simultaneous detection and sizing of small defects in an homogeneous material. However, such a tomographic reconstruction problem is challenging when only limited-angle projections are available, since it is then very hard to recover the geometry of the defects along the average direction of the rays. As far as in situ NDE is considered, this unfavorable observation context is often a direct consequence of the operational constraints imposed on the NDE setup. In [START_REF] Chatellier | Characterization improvement through signal processing application to radiography inspection: 3D recconstruction[END_REF], for instance, the radiographic inspection of in-service pipes in the nuclear industry is considered with a total angular excursion of 30 • , see figure 1.

In limited-angle tomography, faithful reconstructions of the 3D attenuation map cannot be expected without adequate prior constraints on the solution [START_REF] Hanson | Bayesian approach to limited-angle reconstruction in computed tomography[END_REF][START_REF] Rangayyan | Algorithms for limited-view computed tomography: An annotated bibliography and a challenge[END_REF][START_REF] Quinto | Tomographic reconstructions from incomplete data -numerical inversion of the exterior Radon transform[END_REF]. One possible approach is to adopt a parametric model for the shape of the defect (for instance, that it is spherical), and to estimate the few corresponding parameters, say, in the leastsquare sense [START_REF] Rossi | Reconstruction from projection based on detection and estimation of objects -Part I and II: Performance analysis and robustness analysis[END_REF]. However, this approach relies on strong assumptions, which are not natural in the NDE context, where the number and/or the shape of the hypothetical defects are not known. In [START_REF] Benedetti | Multicrack detection in two-dimensionnal structures by means of GA-based strategies[END_REF], a genetic-algorithm strategy is proposed to handle an unknown number of defects, but the constraint of a parameterized shape remains. A potentially more flexible extension is obtained using a description of objects in terms of contours. Again, parametric contour based methods [START_REF] Hanson | Tomographic reconstruction based on flexible geometric models[END_REF][START_REF] Soussen | Polygonal and polyhedral contour reconstruction in computed tomography[END_REF] hardly handle the case of an unknown number of contours, at completly unknown positions. Fully non-parametric modeling of the contour via a level-set formulation [START_REF] Sethian | Level set methods and fast marching methods[END_REF] yields a more suited approach. This leads to gradient-based iterations that aim to compute a local minimizer of the 3D partitioning problem [START_REF] Bruandet | 3D tomographic reconstruction of binary images from cone beam projections: a fast level set approach[END_REF]. A large part of the litterature rather assumes a voxel decomposition for the attenuation map. Then, the tomographic reconstruction x is usually defined as the minimizer of a penalized least-square criterion, such as min

x y -Hx 2 + ψ(x; θ) (1a) subject to x ∈ X ⊆ Ê N (1b)
where θ is a vector of hyperparameters. Within this framework, the solution x is clearly a trade-off: whereas the first term favors the solutions that fit the data, the second one discards the solutions that are not compatible with prior information. A usual choice is to consider ψ as the energy of a Markov random field with positive interactions, in order to favor clusters of voxels in the reconstruction [START_REF] Besag | On the statistical analysis of dirty pictures (with discussion)[END_REF]. More specifically, the well-known Ising model is of special interest here, since the spatial information to encode is binary: each voxel either belongs to an inclusion (of negligible attenuation), or to a sound region of the inspected object (of known attenuation). Only few works have been based on the Ising model to tackle NDE issues. Among the noticeable exceptions, let us mention [START_REF] Retraint | Three-dimensional regularized binary image reconstruction from three two-dimensional projections using a randomized ICM algorithm[END_REF], where it clearly appears that a major numerical obstacle must then be faced, of combinatorial nature, and that no simple, satisfying solution is available. Several authors have rather chosen to relax the binary constraint and to cope with NDE in the framework of continuous-valued Markovian reconstruction [START_REF] Chalmond | Tomographic reconstruction from non-calibrated noisy projections in non-destructive evaluation[END_REF][START_REF] Chatellier | Characterization improvement through signal processing application to radiography inspection: 3D recconstruction[END_REF]. The resulting problem is easier to deal from the numerical viewpoint, since gradient-based optimization algorithms are then available. Unfortunately, the regularizing character of the binary constraint is lost, so that the reconstructed defects appear as streched along some directions [START_REF] Chatellier | Characterization improvement through signal processing application to radiography inspection: 3D recconstruction[END_REF]. Let us remark here that the resort to Markov models strengthens such a tendancy to produce oversized defects. This motivates our choice of an independent, identically distributed (iid) binary model for the voxels: the binary constraint is then the only regularizing ingredient, while the iid assumption corresponds to the simplest possible spatial model. Although the resulting optimization problem is still intractable under an exact form, it lends itself to a simple and satisfactory solution. Its main features are as follows:

i) Firstly, a region-of-interest (ROI) is determined through straightforward computations, such that the maximum posterior mode of voxels outside the ROI is zero, i.e., Pr(

X n = 0 | y) > Pr(X n = 1 | y).
This results in a dramatic reduction in the dimensionality of the problem.

ii) Secondly, we rely on a modified version of the single most likely replacement algorithm (SMLR) introduced in [START_REF] Kormylo | Maximum-likelihood detection and estimation of Bernoulli-Gaussian processes[END_REF] to explore efficiently the binary configurations in the ROI. Whereas the original algorithm considers a "single site" update in the volume, the block most likely replacement (BMLR) proposed here considers "blockwise" updates within cubes of 2 × 2 × 2 voxels. For limited-angle reconstructions, it is empirically observed that such a strategy allows to escape from local minima where the SMLR would be trapped.

Based on a realistic synthetic NDE problem, our numerical procedure is shown to be very efficient, both in terms of computational load and accuracy of reconstructed attenuation map. Moreover, the proposed method does not depend on any hyperparameter, in contrast with many high-resolution tomographic reconstruction methods.

The paper is organized as follows. In Section 2, the modeling of the tomographic measurement in NDE is introduced, and several Bayesian estimators are envisaged to solve the reconstruction problem This section ends with some remarks concerning the computation of these estimators that motivate the design of a computationally efficient heuristic. Section 3 deals with the construction of the ROI that allows a reduction of the voxels involved in the reconstruction process. In Section 4, the problem of binary object reconstruction within a ROI is stated. Standard binary reconstruction algorithms are introduced and the BMLR is derived. In Section 5, the performances of the BMLR are compared to the standard iterated conditional mode algorithm (ICM) and to a penalized approach with convex potentials Gibbs random fields. Finally, a discussion of this work and some conclusions are presented in Section 6.

Binary constraint for NDE

Let S = {1, • • • , N} be an arbitrary enumeration of the N voxels that constitute the unknown three-dimensional attenuation map x. Under usual assumptions, the tomographic measurement along the m-th ray follows the Beer-Lambert relation [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]Section 4.1.1]

: for all m = 1, • • • , M, y m def = -log τ m τ s = h t m x + e m , (2)
where h m ∈ Ê N + stands for the projection process along the mth ray (see for instance [15, Section 7.1]), τ s and τ m are the photon emission rate of the source and the counting rate of detector m, respectively. e m is a noise term that accounts for the statistical fluctuations in the measurement. Stacking all the measurements into a vector y def = (y m ) yields the following observation model

y = Hx + e ( 3
)
where H is a matrix of size M × N, so that h t m is the mth line of H. In this paper, it will be assumed that both the geometry and composition of the inspected object are known, up to the possible presence of defects, of course. In more precise terms, we restrict ourselves to the case where the object would be perfectly known if it was free of any defect. Let r ∈ Ê N + denote the corresponding defect-free attenuation map, and let us introduce a binary map x ∈ {0, 1} N to encode the geometry of the defects: voxels that belong to a defect are assigned the value '1', whereas the others are assigned the value '0'. Given that the defects are assumed to be composed of air, of negligible attenuation, the actual attenuation map x is made of voxels such that

x m = 0 if x m = 1, r m otherwise, for all m = 1, • • • , M.
In compact notations, we have x = rdiag(r) x, so that a new linear observation model can be deduced from (3):

y = H x + e,
where y = Hry, H = Hdiag(r), and e = -e. Thus, the reconstruction problem can be reformulated as a linear inverse problem with binary unknowns. In the rest of the paper, the tilde signs will be omitted for sake of notational simplicity, so that the attenuation map itself will be considered as binary.

In the sequel, e is assumed to be an outcome of a normal random vector E with zero mean and iid components

E ∼ N (0, σ 2 I) ( 4
)
where I is the identity matrix and σ > 0 is the standard deviation of the noise. Provided that the counting rates at each detector are high and of same magnitude, the Gaussian additive observation model is of sufficient accuracy. These assumptions are realistic in NDE and the Gaussian model will be adopted in this paper. For low signal to noise ratio (SNR), however, a Poisson model is more appropriate to describe the random fluctuations in the measurements, see [START_REF] Sauer | A local update strategy for iterative reconstruction from projections[END_REF] for details. Under assumptions (3) and ( 4), the likelihood of the observation vector reads

f Y |X (y | x) ∝ exp - 1 2σ 2 y -Hx 2 (5)
where "∝" stands for "proportional to". Considering the Bayesian framework, a simple iid model for the prior distribution is introduced:

Pr(X = x) = N n=1 Pr(X n = x n ) = p k 1 p N -k 0
where ∀ n ∈ S, Pr(X n = 1) = p 1 , Pr(X n = 0) = p 0 = 1p 1 , and k is the number of voxels that are equal to '1'. Then the posterior probability reads

Pr(X = x | y) = f Y |X (y | x) Pr(X = x) x ′ ∈S f Y |X (y | x ′ ) Pr(X = x ′ ) (6)
where S = {0, 1} N is the set of the binary configurations. The maximum a posteriori (MAP) solution

X MAP def = arg max x∈S Pr(X = x | y) (7)
amounts to minimizing the following quadratic criterion

J (x) def = y -Hx 2 + 2σ 2 µ1 t x (8)
under the binary constraint x ∈ S, where

1 def = (1, • • • , 1) t ∈ Ê N and µ def = log(p 0 /p 1 ).
Beside the MAP solution, the Bayesian setting allows to define alternate estimators. In particular, the posterior mean of voxel n reads

E[X n | Y ] = xn∈{0,1} x n Pr(X n = x n | y) = P n where P n def = Pr(X n = 1 | y).
Following [START_REF] Marroquin | Probabilistic solution of ill-posed problems in computational vision[END_REF], let us also define the marginal posterior mode (MPM) estimator according to:

X MPM n def = arg max xn∈{0,1} Pr(X n = x n | y), (9)
which is easy to deduce from P n :

X MPM n = 1 if P n > 1/2 0 otherwise. ( 10
)
Unfortunately, none of the latter solutions are easily computable. In particular, ( 7) is a binary quadratic programming problem. In the terms and notations of statistical physics, it amounts to minimize the following energy:

H(s) = - (m,n)∈E J mn s m s n + n∈S F n s n , (11)
where

s n = 2x n -1 ∈ {-1, +1} are Ising spins, J mn = -h t m h n , E = (m, n) ∈ S 2 , m < n, h t m h n = 0
is the set of interacting spins, and

F n = -h t n 2y - m∈S h m + 2σ 2 µ
plays the role of the nth component of an exterior magnetic field. In some particular cases, like the one encountered in binary image denoising [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF], solving such a problem can be done in an acceptable number of operations, that is a polynomial function of N. On the contrary, our case falls within the class of NP-hard problems, for which it is very unlikely that a polynomial algorithm exists [START_REF] Barahona | On the computational complexity of Ising spin glass models[END_REF]. The first reason for NP-hardness is the fact that J mn takes negative values (i.e., the interactions are antiferromagnetic) in a nonzero magnetic field. Moreover, the non-planarity of the neighborhood graph E is another source of NP-hardness [START_REF] Istrail | Statistical mechanics, three-dimensionality and NP-completeness: I. Universality of intractability of the partition functions of the Ising model across non-planar lattices[END_REF].

On the other hand, computing X MPM n requires the computation of partition functions involving summations over S or large subsets of S. More specifically, the evaluation of the key quantity P n can be envisaged as follows. Let S i n def = {x ∈ S : x n = i}, so that

P n = x∈S 1 n Pr(X = x | y) = x∈S 1 n Pr(X = x | y) x∈S Pr(X = x | y) = 1 1 + R n (y)
,

where

R n (y) def = x∈S 0 n Pr(X = x | y) x∈S 1 n Pr(X = x | y) . ( 12
)
Again, the computation of the sums entering the expression of R forms NP-hard problems [START_REF] Barahona | On the computational complexity of Ising spin glass models[END_REF][START_REF] Istrail | Statistical mechanics, three-dimensionality and NP-completeness: I. Universality of intractability of the partition functions of the Ising model across non-planar lattices[END_REF]. Despite the NP-hardness of the problem, some contributors resort to Markov Chain Monte-Carlo (MCMC) technics to solve such combinatory problems. In [START_REF] Marroquin | Probabilistic solution of ill-posed problems in computational vision[END_REF], a Gibbs sampler is used to compute the posterior mean estimator for binary denoising problems. In [START_REF] Robert | Binary vascular reconstruction from limited number of cone beam projections[END_REF], the authors resort to simulated annealing in order to minimize a least-square criterion under binary contraint. By nature, such MCMC technics converge in the asymptotic limit. When strong spatial correlations exist in the posterior probability, an efficient sampling is difficult and the resulting algorithms may ask for unrealistic computation time before a useful estimate is available. In practice, we found that this problem is a major obstacle in designing MCMC based techniques for binary reconstruction in limited-angle tomography. Hence, our main concern in this paper is to develop an efficient heuristic that computes a local maximizer of Pr(X = x | y) in a finite number of iterations. For this purpose, our strategy is twofold. Firstly, we aim at reducing the dimensionnality of the problem by an appropriate selection of the voxels that should be discared in the reconstruction process, see Section 3. Then, a local maximizer of Pr(X = x | y) is computed via a deterministic iterative procedure, see Section 4.

Region-of-interest (ROI) selection

As far as large 3D problems are considered, a first step toward a computationally attractive algorithm is to reduce the number of voxels involved in the reconstruction. As shown in the next section, it can be established that X MPM n = 0 for a large amount of voxels, from costless preprocessing operations involving the data.

MPM based ROI in the Gaussian data case

Although P n is not straightforward to calculate, it admits an easily computable upper bound P n that will be fruitful to manipulate.

Let 0 denote the null vector of length N and e n the nth canonical vector of Ê N . Indeed, the following proposition holds in the case of centered white Gaussian observation noise.

Proposition 1 The posterior probability P n admits the following upper bound:

P n ≤ P n def = 1 1 + R n (y) ( 13
)
where

R n (y) def = Pr(X = 0 | y) Pr(X = e n | y) (14) = exp -∆ n , ( 15
)
with ∆ n def = 1 2σ 2 2h t n y -h n 2 -µ. (16)
Proof: See Appendix A. Interestingly, it is clear from [START_REF] Soussen | Polygonal and polyhedral contour reconstruction in computed tomography[END_REF] that

P n ≤ 1/2 =⇒ X MPM n = 0. ( 17
)
The computation of P n plays a pivotal role in our method since it allows a significant reduction of the effective dimension of the reconstruction problem. Let us introduce the set I of voxels for which implication ( 17) is inactive, i.e.,

I = {n ∈ S : P n > 1/2}.
The latter set defines a ROI in the sense that the reconstruction algorithm will be restricted to I, while the other voxels will be assumed to match their MPM estimation, i.e., X MPM n = 0, ∀n ∈ S \ I. Let us remark that

P n ≤ 1/2 ⇐⇒ ∆ n ≤ 0 ⇐⇒ h t n y ≤ 1 2 h n 2 + σ 2 µ ( 18
)
so that a threshold test on the backprojection H t y of the data allows to determine I, where the threshold value is generally not spatially invariant. Moreover, the test does not depend on the noise variance provided that the binary states are considered equally probable, i.e., p 0 = p 1 . A direct extension of Proposition 1 can be established in the general Gaussian case of a noise with mean m and covariance R. A less trivial but richer extension is provided in Subsection 3.3 to encompass some cases of log-concave likelihoods.

Finally, let us note that Proposition 1 admits a formal counterpart, which yields a sufficient condition to ensure X MPM n = 1. The latter condition is hardly met in our situation, because too few voxels are equal to one, but it would be interesting to implement Condition [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF] in more balanced situations. The proof is omitted, since it is almost a paraphrase of that of Proposition 1.

Proposition 2

The posterior probability P n admits the following lower bound:

P n ≥ P n def = 1 1 + R n (y) where R n (y) def = Pr(X = 1 -e n | y) Pr(X = 1 | y) = exp (-∆ n ) , with ∆ n def = 1 2σ 2 2h t n (y -H1) + h n 2 + µ.

Synthetic example

For illustrative purpose, let us compute the threshhold value in a simple, yet realistic context. The inspected object is made of an homogeneous material, with two homogeneous defects of spherical shape. Its discrete representation contains 64 ×64 ×64 voxels, while each of the two defects spreads over 32 voxels. Some schematic illustrations of the measurement context are provided in Fig. 1. The data y consists in seven 2D projections (each one being a 128 × 128 image) that were corrupted by zero mean iid Gaussian perturbations with standard deviation σ = 0.01, which is close to the maximum amplitude of the projection of a single voxel. Both the ROI I and the backprojection were deduced from this data set. The former was computed with the additional assumption that p 1 = p 0 (i.e., µ = 0). Fig. 2(a) and (b) depict the "central" 2D cross-section of the backprojection H t y and of the upper bound P n , respectively.

In this case, 17.3% of the 64×64 voxels of the central cross-section happen to belong to the ROI, and only 1.5% of the 64 × 64 × 64 voxels of the object belong to it, while the ROI incorporates all the 64 voxels of the defects. This example shows that ROI selection can lead to an impressive reduction of the problem dimension. However, it is worth to stress that such a preliminary operation does not raise the main ambiguity concerning the position of the defects, since the voxels that belong to the ROI are organized along the vertical axis, which is the average direction of the rays. It is only aimed at lowering the dimensionality of the reconstruction problem, not at solving it. However, as far as combinatory issues are involved, the size of the vector of unkowns is known as a crucial parameter.

From [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF], it should be clear that the achievable reduction depends on the number of flaws in the volume, the noise in the data, and the number and the angular excursion of the projections. Therefore, it is not easy to assess quantitative results concerning the size of the ROI. However, our experience indicates that increasing the number of projections and/or their angular excursion decreases the size of the ROI. Likewise, a smaller SNR yields a larger ROI.

Generalization to log-concave likelihoods

Proposition 1 can actually be stated in more general terms, provided that the observations y be independent given the object x, and that the likelihood of y m be a log-concave function ‡ of h t m x.

Proposition 3 Let the likelihood function take the following form

f Y |X (y | x) = M m=1 φ m (h t m x), (19)
where φ m are log-concave functions. Then (13)-( 14) still hold.

Proof: See Appendix B.
The main non Gaussian model we have in mind is the Poisson data model, which is the acknowledged model to describe that the measurements are related to a counting process. In the present case of transmission tomography, the corresponding data likelihood reads

Pr(Y = y | X = x) = M m=1 (z m ) ym y m ! exp (z m ) . (20)
In the latter expression, z m = τ s exp (-ξ m ), where τ s is the emission rate of the source and ξ m is related to the attenuation along the ith ray according to

ξ m = h t m x + g m ,
where g = (g m ) stands for a background component that accounts for the Compton scattering effect. Taking the scattering into account (g m = 0) leads to non-logconcave functions. However, when the scattering can be neglected the Poisson likelihood ( 20) is a product of log-concave functions of x and Proposition 3 applies with the corresponding expression for R n (y)

R n (y) = M m=1 exp (τ s (exp (-h mn ) -1) + y m h mn ) . (21)
Here we shall not investigate the Poisson data model further since our main concern is to consider gamma or X ray tomographic measurements with high counting rates, which are well approximated by the simple Gaussian described in Section 2.

Reconstruction using a BMLR approach

In the case of limited-angle tomography, designing a computationaly attractive algorithm that leads to high quality binary reconstructions is a challenge. Here, the proposed approach is twofold. Firstly, the reconstructed volume is reduced to the ROI as defined ‡ By definition, φ is a log-concave function if log φ is a concave function.

in the previous section. Our aim is then to minimize the penalized least-square criterion [START_REF] Benedetti | Multicrack detection in two-dimensionnal structures by means of GA-based strategies[END_REF] with respect to the binary voxels that pertain to the ROI, while the others are forced to zero. For this purpose, we propose a deterministic procedure called block most likely replacement (BMLR). It is an heuristic procedure comparable to ICM. However, it provides significantly better results than ICM, while being still of a moderate computational cost for the considered NDE problems.

ROI constrained formulation

Let us introduce the set I = {x ∈ S : x n = 0 if n ∈ I}. Hereafter, the assumption x ∈ I is taken for granted. It is then natural to consider the problem of maximizing (8) under the constraint x ∈ I. Let us remark here that such a constrained formulation is not necessarily equivalent to the unconstrained one. Actually, it is easy to check that the two problems admit the same solution if X MAP ∈ I only. Unfortunately, it is far from easy to check the latter condition since X MAP is unknown. Nonetheless, maximizing [START_REF] Benedetti | Multicrack detection in two-dimensionnal structures by means of GA-based strategies[END_REF] under the constraint x ∈ I remains coherent and likely to produce a meaningful solution. Moreover, it can be checked by practice that the proposed maximization procedure with or without the ROI provides close reconstruction results -see Section 5.

ICM and SMLR heuristics

>From the sake of low computational complexity, the ICM procedure introduced in [START_REF] Besag | On the statistical analysis of dirty pictures (with discussion)[END_REF] is of particular interest since it amounts to simple successive scalar relaxation of voxels. Let N def = {n 1 , • • • , n |I| } be an arbitrary permutation of the set I. From any initial solution x (1) ∈ I, the kth iteration of ICM can be formally defined by

x (k,1) = x (k) ℓ = 1 • • • |I|, x (k,ℓ+1) = ζ (k,ℓ) if J (ζ (k,ℓ) ) < J (x (k,ℓ) ), x (k,ℓ) otherwise, x (k+1) = x (k,|I|+1) ,
where ζ (k,ℓ) is a proposition vector defined by

ζ (k,ℓ) def = x (k,ℓ) ⊕ e n ℓ , (22)
and ⊕ stands for the componentwise logical exclusive or.

The ICM iterations can be implemented efficiently in a recursive fashion, as follows. Let x and x ′ be two binary N-length vectors that only differ at index n, i.e., x ′ n = x n ±1. Then, it is easy to establish that

J (x ′ ) = J (x) + h t n (h n ∓ 2ε) ± 2σ 2 µ, (23)
where ε = y -Hx. Moreover, we have also

ε ′ = y -Hx ′ = ε ∓ h n (24)
In the implementation of ICM, Equations ( 23) and ( 24) can be used to compare J (ζ (k,ℓ) ) to J (x (k,ℓ) ) and to recursively compute ε (k,ℓ) = y -Hx (k,ℓ) , respectively. By construction, the ICM ensures that the criterion J never increases. Moreover, convergence occurs after a finite number of iterations, as soon as x (k+1) = x (k) . However, an obvious drawback of this heuristic is that the solution is not invariant with respect to the scanning order in I.

The latter limitation of ICM can be made for by a single most likely replacement (SMLR) procedure, as introduced by Kormylo and Mendel [START_REF] Kormylo | Maximum-likelihood detection and estimation of Bernoulli-Gaussian processes[END_REF] in the field of sparse spike train restoration. According to SMLR, each scan yields the swap of only one voxel, the one that produces the largest decrease of the criterion:

x (k+1) = arg max x∈V(x (k) ) J (x), (25)
where V(x (k) ) = {ζ (k,ℓ) : 1 ≤ ℓ ≤ |I|} is the set of all proposition vectors generated by elementary modifications of x (k) . SMLR shares several features of ICM:

• An efficient implementation is obtained thanks to the recursive expressions ( 23) and ( 24).

• The criterion J never increases, and convergence occurs after a finite number of iterations, as soon as x (k+1) = x (k) .

Although very attractive from the computational side, none of these heuristics provide faithfull reconstructions in our context. In particular, we have found that the set V(x (k) ) is not rich enough to prevent ICM and SMLR to get stuck into bad solutions. To gain efficiency, richer sets of propositions should prioritarily allow neighboring voxels to swap their values, so that the position of the hypothetical defects can adjust more easily along the iterations. A possible solution would be to adapt Chi and Mendel's singlespike-shift detector [START_REF] Chi | Improved maximum-likelihood detection and estimation of Bernoulli-Gaussian processes[END_REF]. We found more appropriate to develop an original variation of it, which is developed in the next section.

The BMLR heuristic

Whereas the SMLR only explores the two possible values of a single voxel at each step, the new heuristic explores the 2 8 = 256 configurations of a basic 2 × 2 × 2 cube of voxels in the ROI. The resulting procedure is called block most likely replacement (BMLR).

Formal expression of the heuristic

Let Ω stands for the set of nonempty intersections between the ROI and the 2 × 2 × 2 cubes. Hereafter, such intersections will be called blocks. Most of the blocks contain eight voxels, but some of those located at the boundaries of the ROI may contain a fewer number of voxels.

In what follows, the compact notation x ω is employed for the subvector {x n , n ∈ ω}, where ω is a set of sites, i.e., a subset of {1, . . . , N}. Let us also introduce the notation ζ (k,ω,b) for the proposition vector defined by

ζ (k,ω,b) ω = x (k) ω ⊕ b, ζ (k,ω,b) I\ω = x (k) I\ω ,
where ω ∈ Ω and b ∈ {0, 1} |ω| . Then, the updating equation ( 25) still holds for the BMLR heuristic, with an extended set of proposition vectors:

V(x (k) ) = ζ (k,ω,b) : ω ∈ Ω, b ∈ {0, 1} |ω| .
Again, the criterion J never increases, and convergence occurs after a finite number of iterations, as soon as x (k+1) = x (k) . However, the examination of each block ω requires 2 |ω| evaluations of the criterion, where the bloc-size |ω| varies between 1 and 8. The resulting computational cost can be prohibitive if the implementation issue is not carefully dealt. The next subsection proposes an efficient solution based on Gray codes. Let us remark that the first 2 n elements are obtained by some permutation of the first 2 n natural numbers.

Efficient implementation using

For each block ω, we propose to explore the 2 |ω| binary configurations in the Gray code order rather than in the usual one. Then, two successive trials b and b ′ respectively correspond to two vectors ζ (k,ω,b) and ζ (k,ω,b ′ ) that only differ at a single index. As a consequence, the recursive expressions ( 23) and ( 24) can be applied to explore the 2 |ω| values of J (ζ (k,ω,b) ) at an affordable numerical cost, for all ω ∈ Ω.

An even more efficient version consists in introducing f n = h t n ε and g mn = h t m h n for all m, n ∈ ω, so that ( 23) can be replaced by

J (x ′ ) = J (x) + g nn ∓ 2f n ± 2σ 2 µ,
while the |ω| scalars f n can be recursively computed using

f ′ m = h t n ε ′ = f m ∓ g mn .
While each of the 2 |ω| applications of ( 23) involves a scalar product, none is required anymore in the new recursion. Only |ω| (|ω| + 1)/2 scalar products are required to compute g mn for all m, n ∈ ω. Table 1 depicts the pseudo-code for the central part of the BMLR scheme. For the sake of brevity, the set of blocks Ω is assumed available from previous calculation, as well as the projection matrix H. For large size problems, H cannot be stored and its entries must rather be recalculated when required. In Table 1, the initial configuration is the null object x 0 = 0. However, a different initialization point could be easily considered. In particular, a faster but less optimal method could provide an initial solution, e.g., ICM, that would be refined using BMLR. For the sake of computational time, such a two-step approach would be of particular interest to reconstruct objects of larger extent than localized defects.

Test on synthetic data sets

The 3D synthetic NDE problem presented in Section 3.2 is now considered in order to test the proposed reconstruction algorithm. The aim is to illustrate the detection, sizing and separation capabilities of the method on the ground of simulations.

The simulation context

The considered reconstruction problem involves two superimposed defects of spherical shape, in an homogeneous material. The volume of interest contains 64 × 64 × 64 voxels, while each of the two defects spreads over 32 voxels. Two different values have been considered for the distance between the two defects: L ∈ {26mm, 12mm} (see Figure 3(a)). The data set y consists in seven projections that were corrupted by an uncorrelated Gaussian noise with standard deviation σ ∈ {0, 0.005, 0.01} (see Figure 3(b)). The geometric dimensions of the problem are given by Figure 1. They correspond to a realistic situation of in-service pipe inspection in the nuclear industry [1, Fig. 1].

Non-binary reconstruction results

Here, a non-binary penalized approach tested in [START_REF] Chatellier | Characterization improvement through signal processing application to radiography inspection: 3D recconstruction[END_REF] is used as a comparative tool. It consists of considering an optimization problem of the form (1) with X = Ê N + and

ψ(x; θ) = λ i∼j φ(x i -x j ; δ) + α N n=1 |x n | , (26)
where λ, α ≥ 0, i ∼ j stands for all pairs of neighboring sites in the first-order neighborood system (hence, each voxel has six neighbors) and φ(t; δ) is a Huber function [START_REF] Huber | Robust Statistics[END_REF]. The role of such a penalty function is to favor piecewise homogeneous solutions.

In practice, the penalized approach requires the tuning of the three hyperparameters θ = (λ, δ, α) t . Here, they have been empirically adjusted according to qualitative appreciation. The corresponding reconstruction results are depicted by Fig. 4. As expected, the estimated attenuations are underestimated, while the extension of the defects is overestimated along the vertical axis. In particular, the result does not allow to clearly discriminate the two spheres in the case L = 12mm. The penalized least-square criterion being convex, the global minimizer can be computed from an initial guess x (0) by an appropriate iterative descent algorithm. Here, x (0) has been chosen as the backprojection of the data and the minimization step was performed using a projected-gradient method, i.e., a straightforward adaptation of steepest-descent that handles the positivity constraint in a natural way [25, p.203]. However, this algorithm inherits the slow asymptotic convergence of a gradient method. Table 2 shows that several hours are required to reach the adopted convergence criterion ψ(x (k-1) ) -ψ(x (k) ) < 10 -6 , where k ∈ AE the iteration number. Let us remark here that the same constrained minimization problem could probably be adressed using more efficient, but more complex algorithms (see for instance [START_REF] Johnson | A primal-dual method for large-scale image reconstruction in emission tomography[END_REF]).

Binary reconstruction results

In this section, the proposed BMLR algorithm is compared to the ICM. In all cases, the initial state is chosen as a defect-free object, i.e., x (0) = 0. Both heuristics require the tuning of the hyperparameter µ = log(p 0 /p 1 ). Since the defects are expected to constitute a very small part of the total volume, a strict Bayesian viewpoint would lead to choose prior probabilities p 0 and p 1 = 1p 0 such that p 0 ≫ p 1 . For two distinct reasons, we rather set p 0 = p 1 .

• On the one hand, false positives are preferable to false negatives, which justifies to overvalue p 1 /p 0 . This could be reformulated more rigorously according to the Bayesian cost theory, by replacing the MAP estimator by another Bayesian estimator that accounts for distinct costs for false negatives and false positives, respectively c 0 and c 1 . Then our choice corresponds to the assumption that p 0 /c 0 = p 1 /c 1 .

• On the other hand, choosing p 0 = p 1 leads to a fully unsupervised version of BMLR. The reconstruction results are shown in Figs. 5 and6 for the ICM and the BMLR algorithms, respectively. Even in the noiseless case, the ICM fails to provide a faithful reconstruction in the case of two close spheres. In comparison, the BMLR algorithm achieves a perfect reconstruction in all cases of high or moderate SNR.

On the other hand, the number of false positive voxels is subject to a threshold effect: when the noise standard deviation reaches the magnitude of the projection of a single voxel (whose maximum value is σ th = 1/64 ≈ 0.015 in the present simulation), many spurious, isolated voxels appear in the reconstruction. This is clearly visible on Fig. 6 in the case σ = 0.01. In this case, two distinct strategies may be adopted to reduce the number of false positive voxels. The first one rests upon the fact that the BMLR reconstruction remains acceptable, let alone those spurious voxels. As far as real defects can be assumed to extend over several neighboring voxels, the idea is then to remove isolated voxels by simple postprocessing. For sake of computational burden, it is interesting to also discard isolated voxels from the ROI. In the considered cases where σ = 0.01, a large amount of iterations are spent to activate voxels that are isolated in the ROI, and can thus be saved. The resulting reconstruction is depicted in the first column of Fig. 8.

An alternate strategy is to introduce a penalization of active voxels by assigning a positive value to the regularization parameter λ def = 2σ 2 µ in criterion [START_REF] Benedetti | Multicrack detection in two-dimensionnal structures by means of GA-based strategies[END_REF]. The tuning of λ has a great impact on the final quality of the 3D attenuation map. As illustrated in Fig. 9, a correct tuning of λ significantly reduces the number of spurious voxels and produces a satisfactory reconstruction of the flaws. However, over estimating λ leads to false negative voxels, so that the tuning of λ must be cautiously handled. The number of iterations and the CPU time required for convergence for the ICM and BMLR heuristics are shown in Table 2. Not surprinsingly, ICM is a very fast method. However, BMLR remains very attractive in terms of computing time. The step of building the ROI is not expensive since it results from a simple thresholding of the backprojection map, see [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF]. Depending on the SNR and on the positions of the spheres, the ROI varies from 1.1% to 1.6% of the total number of voxels in the considered examples. The total number of iterations is roughly proportional to the number of activated voxels in the solution. Obviously, the latter depends on the volume of the defects, which is expected to be small in NDE. However, it is also dependent on the SNR: the BMLR heuristic converges in a dozen of iterations if σ < σ th whereas a few hundreds are required if σ > σ th . In the former case, the convergence is fast and the solution is reached within two seconds. In the latter case, most of the iterations activate only one isolated site of the ROI and the time needed for convergence increases significantly. As shown in Table 2, removing isolated sites from the ROI allows to speed- up the convergence of the BMLR in this case. As stated earlier in this section, this pre-processing of the ROI is natural provided that flaws extend over more than a single voxel in the volume.

In principle, the BMLR heuristic can provide different solutions whether the ROI is considered or not. In order to test the variability of the solution, we also ran the BMLR heuristic with I = S, i.e., we ignored the ROI. For high to moderate SNR values, the results were strictly identical to the previous ones, while only marginal differences appeared for the low SNR case. Finally, let us remark that the BMLR heuristic asks for an unrealistic computation time of several hours when the ROI is ignored, or if all the voxels belong to the ROI. In such a defavorable situation, the BMLR could still be efficiently implemented using a parallel structure, since the blocks can be explored independently one from each other.

BMLR: robustness with respect to modeling errors

In the previous subsection, the data are noisy projections of discretized objects. Obviously, real objects are continuous and (3) only constitutes an approximate model. Therefore, we found appropriate to further investigate the robustness of the method by considering projections from perfectly spherical objects, given that closed-form expressions are available for the projection of ellipsoids. The reconstructions provided by the BMLR are shown in Figure 7. In the case of distant objects, a good reconstruction is obtained whatever the SNR -up to the fact that spurious, isolated voxels must be removed for σ = 0.01, as shown on the second column of Figure 8. In the case of closer objects, the quality of the reconstruction deteriorates slightly, but the BMLR still provides meaningful results, since the two objects are still clearly distinguishable one from each other. does not efficiently eliminate spurious voxels and over regularization (right) cancels most of the voxels in the flaws, a correct tuning of λ (center) produces a satisfactory reconstruction.

Conclusion

Spatially interacting models such as Markov random fields are usually employed within a Bayesian approach to reconstruction. Here, we rather resorted to a simpler iid model, and we proposed an efficient binary reconstruction algorithm, both in terms of robustness and computational burden. For high to moderate SNR, the binary constraint is sufficient to regularize the inversion of the tomographic problem. For low SNR, however, the number of false positive voxels increases dramatically. A simple strategy to deal with this problem is to cancel isolated voxels. Another strategy is to perform a penalization of active voxels by assigning a positive value to the regularization parameter λ = 2σ 2 µ in criterion [START_REF] Benedetti | Multicrack detection in two-dimensionnal structures by means of GA-based strategies[END_REF]. However, the supervised tuning of λ must then be handled. With a view to reconstruction from real data, several comments can be made about the physical model considered here. Firstly, the additive Gaussian noise assumption may not be appropriate for low counting rates τ m . In such cases, the Poisson model [START_REF] Istrail | Statistical mechanics, three-dimensionality and NP-completeness: I. Universality of intractability of the partition functions of the Ising model across non-planar lattices[END_REF] provides a more accurate description of the data. The ROI should then be computed using [START_REF] Robert | Binary vascular reconstruction from limited number of cone beam projections[END_REF], and straightforward modifications should be brought to the BMLR algorithm.

Secondly, it should be noted that the Beer-Lambert relation ( 2) is only valid for a monochromatic emission of photons under a non-diffractive propagation assumption. In particular, Compton scattering is then negligible, which is not necessarily true for high energy gamma-ray sources (beyond 100 keV). In the latter case, further investigations would be required to define an appropriate ROI, for instance on the basis of the scattering model developped in [START_REF] Sauer | Bayesian estimation of 3-D objects from few radiographs[END_REF].

FilmFigure 1 .

 1 Figure 1. Left -3D representation of the instrumental setup. Right -Schematic representation of the measurement context for the synthetic problem considered in this paper. Note that all the dimensions are drawn from a real NDE setup in the nuclear industry, cf. [1, Fig. 1].

Figure 2 .

 2 Figure 2. Central 2D section of the backprojection (left) and of the map of the upper bound P n (right), repectively, obtained from noisy projections (σ = 0.01), for a pair of superimposed defects (L = 26mm).

  Gray codes A Gray code (or reflected code) is a binary numeral system where two successive values differ in only one digit [23, Section 20.2]: (0, 1, 11, 10, 110, 111, 101, 100, 11000, 11001, . . .)

Figure 3 .

 3 Figure 3. Left -Pairs of voxelized spheres considered in the data generation process, in configuration separated (L = 26mm) and close (L = 12mm), respectively. Right -Corresponding noisy projection from source S6 with noise standard deviation σ = 0.01.

Figure 4 .

 4 Figure 4. Reconstructions from noisy projections of a pair of voxelized spheres using a penalized approach.

Figure 5 .

 5 Figure 5. Reconstructions from noisy projections of a pair of voxelized spheres using an ICM approach.

Figure 6 .

 6 Figure 6. BMLR reconstructions from noisy projections of a pair of discretized defects.

Figure 7 .

 7 Figure 7. BMLR reconstructions from noisy projections of a pair of continuous defects.

Figure 8 .

 8 Figure 8. For the lowest SNR values, removing isolated voxels in the final BMLR reconstruction provides a better solution to both discretized (left) and continuous (right) defect problems.

λ = 1 . 7 × 10 - 4 λ = 1 . 7 × 10 - 3 λ = 1 . 7 × 10 Figure 9 .

 171041710317109 Figure 9. Impact of the regularization parameter λ = 2σ 2 log(p 0 /p 1 ) on BMLR reconstruction at the lowest SNR (i.e., σ = 0.01). Whereas under regularization (left) does not efficiently eliminate spurious voxels and over regularization (right) cancels most of the voxels in the flaws, a correct tuning of λ (center) produces a satisfactory reconstruction.

Finally

  , term-to-term summation of inequalities (B.1) and (B.2) allows to conclude that φ m (a) φ m (d) ≤ φ m (b) φ m (c).

Table 2 .

 2 Number of iterations and CPU time for the penalized approach and the ICM and BMLR heuristics on a Pentium III 2.4 Ghz (586 series) with 1Go RAM.

	Algorithm	L [mm]	σ	# it. CPU time
		26	0	333 3h 5min
			5.10 -3	97 1h 25min
	Penalized approach	12	10 -2 0	198 2h 55min 325 3h 42min
			5.10 -3	356 5h 11min
			10 -2	410 5h 58min
		26	0	9 0.408s
			5.10 -3	11 0.452s
	ICM	12	10 -2 0	11 0.501s 7 0.306s
			5.10 -3	8 0.331s
			10 -2	8 0.333s
		26	0	30 1.57s
			5.10 -3	27 1.55s
	BMLR	12	10 -2 0	346 16.94s 19 1.18s
			5.10 -3	19 1.23s
			10 -2	339 17.08s
	BMLR	26	0	30 1.57s
	(with		5.10 -3	27 1.53s
	isolated voxels discarded from the ROI)	12	10 -2 0 5.10 -3 10 -2	173 8.37s 19 1.17s 19 1.29s 170 8.43s

Acknowledgments

This work was partially supported by Electricité de France (EDF) under contract number E961/P10AY0/E3.

end for

which readily gives expression [START_REF] Sauer | A local update strategy for iterative reconstruction from projections[END_REF]. Now let us show that

for all couples (x, x ′ ) ∈ S 0 n × S 1 n such that x ′ = x + e n : given ( 5), the latter inequality holds, since

Given [START_REF] Bruandet | 3D tomographic reconstruction of binary images from cone beam projections: a fast level set approach[END_REF], inequalities (A.1) jointly imply that (13) holds.

Appendix B. Proof of Proposition 3

In the more general setting of Proposition 3, let us show that (A.1) is still valid for all couples (x, x ′ ) ∈ S 0 n × S 1 n such that x ′ = x + e n . Given [START_REF] Barahona | On the computational complexity of Ising spin glass models[END_REF], this amounts to prove that


S Leclerc

M Petryk

D Canet

J Fraissard
email: jacques.fraissard@upmc.fr

Competitive Diffusion of Gases in a Zeolite Using Proton NMR and a Slice Selection Procedure

Keywords: co-diffusion in zeolite, co-adsorption in zeolite, benzene diffusion in ZSM5, hexane diffusion in ZSM5

The study of the co-diffusion of several gases through a microporous solid and the resulting instantaneous distribution (out of equilibrium) of the adsorbed phases is particularly important in many fields, such as gas separation, heterogeneous catalysis, etc. Classical 1 H NMR imaging is a good technique for the visualization of these processes but, since the signal obtained is not specific for each gas, this requires that each experiment be performed several times under identical conditions, and each time with only one incompletely deuteriated gas. In contrast, we have proposed a new NMR imaging technique (based on the so-called NMR lift) which gives a signal characteristic of each adsorbed gas. It can therefore provide directly, at every moment and at every level of the crystallite bed, the distribution of several gases competing in diffusion and adsorption.

Introduction

Knowledge of the co-diffusion and co-adsorption coefficients of reactants and products is essential when a heterogeneous catalysis reaction is performed by means of a gas flowing through a porous catalyst bed. It can be even more important in a fluidized bed where the contact times between the catalyst pellets and the reactants are very short, which means that one hardly ever knows whether the pellet is used completely or only superficially. Under these experimental conditions the system is never at adsorption equilibrium. The distribution of the various reactants adsorbed on the catalyst is very inhomogeneous and, moreover, very variable from one reactant to another. However, calculations on the kinetics are generally performed by assuming the ideal case of a homogeneous assembly, both for the totality of the catalyst and for the distribution of all the reactants, which makes these calculations, even when they are mathematically rigorous, rather approximate as regards reality. Let us recall, for example, that in the case of a reaction following the Langmuir-Hinshelwood mechanism, the reaction rate is given by the product of the local concentrations of all the reactants, integrated over the entire catalyst. It is therefore necessary to be able to determine at every moment the diffusion coefficient of each reactant in the presence of the others and its instantaneous distribution along the length of the catalyst bed. Classical 1 H NMR imaging is a good technique for the visualization of these processes but, since the signal obtained is not specific for each gas, this requires that each experiment be performed several times under identical conditions, and each time with only one incompletely deuteriated gas. This is very time-consuming and becomes very expensive when completely deuteriated molecules are required. Despite these difficulties it has been previously demonstrated that the distribution of the adsorption of benzene and hexane flowing through a ZSM5 bed is very heterogeneous [1]. To remedy the drawbacks of classical imaging we have proposed [2] a new NMR imaging technique which can be used for following the diffusion and adsorption of a gas in a microporous bed [3]. The sample is displaced vertically, step by step, relative to the detector during the adsorption of the gas; the detector is a very thin coil. The bed is assumed to consist of n very thin layers of solid, and the probed region is limited to each layer, so that the variation of the concentration of gas absorbed at the level of each layer is obtained as a function of time. This technique permits the determination of all the parameters of the system at every moment, at every position in the solid and in each crystallite. The complex calculations used to treat the diffusion of benzene alone through a ZSM5 bed have already been published [3,4]. An interesting feature of this technique is its ability to visualize directly the co-diffusion of several gases. Indeed, the NMR signal is quantitatively characteristic of the adsorbed gas. It can therefore provide directly, at every moment and at every level of the crystallite bed, the distribution of several gases competing in diffusion and adsorption. We only present here the interest of this technique, without going into the calculation of the resulting diffusion and adsorption parameters.

Experimental

Our aim is to determine the evolution of the NMR spectrum along the length of a heterogeneous sample in the z direction of the magnetic field 0 B . This original method of imaging has been presented in [1]. The idea is to select successively several sections at different positions of the sample along the z axis. For this purpose a mechanical device (NMR lift) was built to displace the sample tube (external diameter 5 mm) vertically step-bystep inside the NMR probe. In this way one obtains for each position a spectrum corresponding to the composition of the sample at this point. The experiment presented here concerns, as a first example, the co-diffusion of benzene and hexane through a ZSM5 bed (length 15 mm). Figure 1 shows the sample-holder tube which is moved vertically, up or down, opposite the very thin detector. The homogeneous liquid solution consists of two equal volumes of benzene and hexane. It is in equilibrium with its gas phase at 25°C. The two gases begin to diffuse in the zeolite when the glass partition is broken. It is clear that in such an experiment we cannot control all the physical parameters. In order to minimize the temperature variation when the glass wall is broken and the gas comes into contact with the zeolite, the volume under vacuum above the solid is very small and practically negligible compared to the much greater volume of the upper ampoule containing the liquid or gaseous benzene. Moreover, the area of liquid is relatively high, and the diameter of the glass wall to be broken is not small for not disturbing the diffusion of the gas. Finally, because of problems of pollution and of the construction of the NMR probe, we had to use a closed sample-holder tube and not a completely open one. This restriction in no way diminishes the interest of the method.

In accordance with the principle mentioned above, the zeolite sample is displaced vertically in the magnet during the adsorption by 0.5 mm steps. The signal of the adsorbed gas is detected for each position, rapidly so as to obtain the best possible precision. When the sample reaches the uppermost or the lowest position it starts again in the opposite direction, and continues in this way until adsorption equilibrium is reached. This procedure thus provides the amount of gas adsorption as a function of time for each level in the adsorbent. This makes it possible to determine a diffusion coefficient at each point of the solid [2]. This is essential for tow studies that we are planning : (i) the adsorption of gas in beds consisting of layers of different solids, and (ii) the competitive diffusion of several gases. The latter study will be facilitated by the fact that one detects directly the real characteristic NMR signals. Right: scheme of the narrow which is zone monitored by the detector.

Results

Figure 2 shows the 1 H NMR spectra of benzene and hexane at two positions, at 6 and 14 mm (near the bottom of the tube) and after two adsorption times, 34 and 280 minutes. After 34 min. there is still a high concentration gradient the length of the sample. The equalized intensities of the two spectra after 280 min. show that the system has reached equilibrium. To compare the concentrations of adsorbed benzene and hexane molecules, a coefficient, 6/14, has to be applied to the spectrum of the latter to take into account the numbers of H atoms in each compound. This is displayed in Figures 3 and4 Figure 3 presents the benzene and hexane concentrations along the sample, from a height of 6 mm down to 14 mm, for each diffusion time. We have not considered neither the first diffusion time nor the first layers very close to the gas phase, because it takes 5 to 6 min. from the moment when the glass is broken (therefore the beginning of diffusion) to get the first spectra. It is obvious that after such a time the first layers are saturated, but this restriction does not affect the calculation of the diffusion parameters (not performed here). This figure displays the time dependence of the concentration gradient of each of the gases adsorbed along the catalyst. The signal intensity scale is the same for both components. Figure 4 compares the evolution, as a function of time, of the benzene and hexane concentrations, at different levels of the sample. It reveals particularly well, under the chosen experimental conditions, the hindrance of benzene upon the diffusion of hexane, and this at every moment.

Conclusion

The choice of the two diffusing gases, benzene and hexane, is not what is most important. In the actual experiment, the system reaches equilibrium and adsorbat homogeneity along the catalyst after a certain time because we had to use a tube closed at the bottom to prevent external pollution. The system will not reach neither equilibrium nor homogeneous distribution when the gases go through a reactor open at both ends, as is generally the case for a reaction. The main result is the possibility, for the first time in a single experiment, of probing at every moment the concomitant distribution of several gases co-diffusing in a porous solid.

Fig. 1 .

 1 Fig. 1. Left: sample-holder bulb containing the liquid phase in equilibrium with the gas phase;Right: scheme of the narrow which is zone monitored by the detector.

  Figure2shows the 1 H NMR spectra of benzene and hexane at two positions, at 6 and 14 mm (near the bottom of the tube) and after two adsorption times, 34 and 280 minutes. After 34 min. there is still a high concentration gradient the length of the sample. The equalized intensities of the two spectra after 280 min. show that the system has reached equilibrium. To compare the concentrations of adsorbed benzene and hexane molecules, a coefficient, 6/14, has to be applied to the spectrum of the latter to take into account the numbers of H atoms in each compound. This is displayed in Figures3 and 4.

Fig. 2 .

 2 Fig. 2. 1 H NMR spectra of benzene (B, 6 atoms) and hexane (H, 14 atoms) at 6 and 14 mm (from the top of the tube) after 34 and 280 minutes.(The numbers on the abscissa axis are only indicative).

Fig. 3 .

 3 Fig. 3. Benzene and hexane concentrations along the sample, from 6 mm to 14 mm, for each diffusion time.

Fig. 4 .

 4 Fig. 4. Evolution, as a function of time, of the benzene and hexane concentrations, at different levels of the sample.

Figure 1 .

 1 Figure 1. Left: Sample-holder bulb containing the liquid phase in equilibrium with the gas phase; right: Schema of the narrow zone monitored.

Figure 2 . 1 HFigure 3 .

 213 Figure 2. 1 H NMR spectra of benzene (B, 6 atoms hydrogen) and hexane (H, 14 atoms hydrogen) at 6 and 14 mm (from the top of the tube) after 34 and 280 minutes.(The numbers on the abscissa do not correspond exactly to the chemical shifts)

Figure 4 .

 4 Figure 4. Evolution, as a function of time, of the benzene and hexane concentrations, at different levels of the sample.