Multiexponential relaxation and local diffusion cells distribution

Uniform PENalty inversion algorithm (UPEN)














Multiexponential relaxation 
and local diffusion cells distribution


In bulk water the relaxation of the 1H nuclei magnetization is well described by a single exponential function. For water in high S/V systems the relaxation is enhanced, and shorter relaxation times are observed. 
The semi-logarithmic plot
s of the decay of the nuclear magnetization in two water-saturated model systems, each one characterized by a narrow distribution of  pore sizes,  are straight lines, with faster decay for the sample with the smaller pores. Both systems exhibit faster decay than for bulk water.

The basic phenomenon used to interpret this effect is the diffusion of a water molecule from the centre of a pore to the walls, where the relaxation is enhanced. The strength of the surface effect is usually represented by the parameter ρ (having the dimensions of a velocity). If diffusion is sufficient to maintain the nuclear magnetization nearly constant over the pore space, including at the surfaces, and if bulk relaxation is negligible (conditions usually satisfied in porous rocks), the relaxation time is inversely proportional to S/V and the simple relationship 1/T1,2 = ρ S/V  will apply. Here the subscripts 1 and 2 of the relaxation time stand for the longitudinal (T1) and transverse (T2) components of the nuclear magnetization, respectively, of 1H nuclei.

What will happen in real permeable porous media, made of networks of interconnected pores, of different shapes, and with distributions of characteristic lengths?

Relaxation represent an average over regions that can be much larger than an “individual pore”. If  diffusion is fast enough to maintain magnetization nearly constant within distances of diffusion in a local relaxation time, that is, within the local diffusion cell, the local relaxation rate is proportional to an  S/V averaged over this local volume. It is often useful to regard the local V/S as an operational definition of “pore size”. 


If the pore-space is heterogeneous in ρ and/or S/V only on a scale shorter than local diffusion lengths, which implies uniform local diffusion lenghts and local relaxation times, the relaxation of the sample will be single-exponential. Otherwise, the distribution of local S/V will be reflected by a multiexponential relaxation [J. Appl. Physics 1996, 3656-3664].

The model of local relaxation averaging over local diffusion
cells was well validated [Magn. Reson. Imaging, 12, 1994, 191-195] by an experiment of progressive desaturation of an initially fully-saturated rock sample.  If we saturate a sandstone sample with water, we get a distribution covering about three decades, from a ms to one second. After progressive desaturation by centrifugation of the sample we have a clear visualization of the desaturation process: most of the water leaves the larger pores, so we observe a progressive vanishing of the amplitudes at high T1 values. The  progressive shift of some of the signal to shorter T1 values reflects less mixing by diffusion between adjacent pores of different sizes, one or both of which is partly emptied. Furthermore the larger pores have probably lost a larger fraction of fluid. Thus, both small and large pores may have relaxation time decreased.





























Inversion from signal time to relaxation times


The problem is then to invert the multiexponential function in order to obtain the distribution of relaxation times.

The function to minimize is the squared error of fit plus a penalty term. Usually some smoothing is implemented in the distributions by means of a penalty  function in order to avoid excess variation. We have proposed an algorithm (Uniform-PENalty, UPEN)) to obtain distributions of relaxation times, tested on rocks and other porous materials, as well as artificial data, where we could see distributions more than three decades wide. UPEN uses a smoothing coefficient varying with the relaxation time and determined by iterative feedback in such a way that the smoothing penalty, rather then the coefficient, is roughly uniform.

UPEN: Uniform PENalty inversion of multiexponential decay data



Minimize  (error)2 + penalty

The first line is weighted (error of fit)2; the second is the penalty.
Bi is the weighting factor, 
gk the computed amplitude,
ti is data time, 
Tk is relaxation time, 
si is signal amplitude. 
Ak is amplitude penalty coefficient (=0 for times covered by data),
Ck is curvature penalty coefficient (with negative feedback applied to the penalty through ck).






ßp is the slope (pendenza) compliance parameter.
ßc is the curvature compliance parameter.
ßa is the amplitude compliance (for outside data range).
ß0 is the compliance floor (not critical, but not zero).
ß00 (normally = 1) multiplies the other ß 's.
Q is data point time spacing in Nepers.

UPEN tends to avoid more maxima and minima in the distribution than demanded by the data, minimizing the appearance of separate populations in the computed distributions to the extent permitted by the data.

This is an example of a carbonate rock saturated by water. UPEN gives a sharp line that is not broadened more than is consistent with the noise and in the same distribution shows a tail decades long without breaking it up into several peaks not required by the data to be separate. On the contrary, other solutions obtained by fixed smoothing coefficient or by singular-value decomposition (with substantially the same rms error of the fit) in any case tend to broaden the major peak  and break the tail into more peaks, which might be interpreted as physically meaningful resolved compartments.











Errors of Fit, Noise, Quality of Fit



The Ei 's are the errors of fit.  For T2 data, signal amplitude si may represent a window average of Bi echoes.  Rr (R for noise: Rumore, Rauschen) is an estimate of individual-echo error of fit, using a compromise weighting to give equal weights at long and short times, rather than Σ(BiEi)2/ΣBi to give the best statistical estimate of the ensemble average echo noise R.  This compromise weighting is used for all the R's.

If there is a slowly-varying error of fit in addition to random noise, one may get a better estimate of the noise from Rv (v for variation), computed from differences of the E 's for second-nearest neighbors, since a non-random slowly-varying error of fit would not greatly affect these differences.

If there are random errors only, Rr should be somewhat less than Rv , because the fitting minimizes the errors, rather than differences of errors.  Having separate estimates of the random noise and the overall error of fit can warn us if the data are not, or cannot be, adequately fit.  Thus, the quality-of-fit parameter

Rrv = ln ( Rr / Rv )

should normally be negative, as it usually is for artificial data corresponding to positive exponential components plus random noise.  A value of Rv=0.05 usually indicates data problems, and values over 0.1 normally indicate serious data problems.














Resolution of Two Equal Lines


E2 is the least maximum absolute error of fit to a rectangular distribution of relaxation times by 2 sharp lines, a factor of Y apart.  The data time spacing in Nepers is Dq, which for T2 data without discarded points is TEE/T, where T is relaxation time of the minimum between the peaks, and TEE is the spacing of the echoes employed.  If the peaks are of comparable size and are the only significant features on the distribution, SNR is the signal-to-noise ratio.  If the peaks are disparate, only the signal in the smaller peak is counted.  Additional major features make resolution more difficult.  The above SNR value is for marginal resolution; several times higher SNR is required for firm resolution.  If the two tentative peaks do not have zero inherent line width, Y must be reduced by about the mean halfwidth (on a logarithmic time scale).












Linewidth due to the Noise



where w is the halfwidth of the line and SNR counts only the signal in the line.  Nearby features widen the line.









Significance Criterion for Improvement in Fit



Under some conditions, to suggest a significant change in the quality of a fit, the relative change in error of fit should be at least


with the same adaptation to T2 data as mentioned above.  N is the number of (possibly windowed) data points used in the computation.