Measuring adsorption, diffusion and flow in chemical engineering: applications of magnetic resonance to porous media

Magnetic resonance (MR) techniques are increasingly used to improve our understanding of the multi-component, multi-phase processes encountered in chemical engineering. This review brings together many of the MR techniques used, and often developed specifically, to study chemical engineering systems and, in particular, processes occurring within porous media. Pulse sequences for relaxometry, pulsed field gradient measurements of diffusion, imaging and velocimetry measurements are described. Recent applications of these MR pulse sequences to microporous, mesoporous and macroporous structures are then reviewed. Considering the microporous and mesoporous systems, we focus attention on studies of rock cores, manufactured materials such as cement and gypsum plaster, and catalysts. When considering macroporous structures, the transport through packed beds of particles typical of fixed-bed catalytic reactors is reviewed; a brief overview of the increasing research interest in gas–solid fluidized beds is also presented. We highlight the field of sparse k-space sampling as an area that is in its infancy and suggest that, combined with Bayesian methods, it will offer new opportunities in both extending the application of high-field MR techniques to chemical engineering and increasing the range of measurements that can be carried out using low-field hardware.


Introduction
Since the early days of applying magnetic resonance (MR) to problems in applied physics and chemistry (Packer 1967, Derouane et al 1972, Horsfield et al 1989, Webb et al 1989, Smith and Strange 1996, Watson and Chang 1997, Miller 1998, MR techniques have found increasing application in mainstream chemical engineering (Gladden 1994, Gladden 2003, White and Truitt 2007, Hunger 2008. Many of these more recent applications have required the development and implementation of novel pulse sequences to address the challenges associated with magnetically heterogeneous, short relaxation time systems. This review summarizes the technical developments that have taken place mostly in the last 10 years with a particular focus on the applications to porous media, and illustrates the type of data that can be obtained now. We focus on methods developed by the Cambridge group but put these in the context of the wider research activity in this field. Developments have occurred in fast data acquisition across the whole range of MR techniques we discuss. In parallel with this is the need to implement MR to discriminate between different phases of matter and between different chemical species. The review of techniques available to us for exploring porous media in chemical engineering is divided into three sections: (i) relaxation time analysis of liquids in porous media, (ii) flow propagator measurements and (iii) magnetic resonance imaging and velocimetry (also referred to as flow imaging or flow mapping). In practice, we often integrate aspects of these different measurements into a single acquisition. We split the review of the applications of these pulse sequences into two broad themes. The first theme considers microporous and mesoporous systems with a particular focus on studies of rock cores and model porous media, materials such as cements and gypsum, and catalysts. Relaxometry and propagator measurements are usually the techniques of choice for studying these systems; however, the precise information required from each study can be quite different, depending on the process application. The second theme is that of macroporous systems, in which we focus on fixed-bed reactors and gas-solid fluidized beds. Here, magnetic resonance imaging (MRI) and flow mapping techniques are used more commonly. In taking these methods forward to ever shorter data acquisition times, we note the emerging area of sparse sampling in 3 applications to chemical engineering systems and the implementation of Bayesian methodology in data acquisition strategies. The opportunities for undertaking integrated experimental and numerical simulation studies of transport in porous media are also highlighted, not only in the context of validating the numerical codes but also in using the codes to explore the physical phenomena determining the adsorption, transport and reaction processes occurring within the porous medium of interest.

Technique development
The main drivers in technique development have been the desire to acquire quantitative data, often with the additional requirement that the system of interest is in an unsteady state. Therefore, we require the data acquisition to be as fast as possible. It is also important to consider what we mean by 'quantitative'. In the context of transport (diffusion and flow) and reactions in porous materials, we expect to get good agreement-within a defined error bar-on voidage (porosity), total flow rates and chemical compositions determined by traditional porosimetry, gravimetric analysis and gas chromatography measurements, respectively. Considering the concept of 'fast' measurements, MRI experiments that gather the full k-space raster have now probably reached the limit of these data acquisition strategies. Interest has therefore shifted to a focus on sparse sampling of k-space, and this is an ever growing area of research. In particular, we will highlight the use of compressed sensing (CS) methodology and Bayesian analysis of experimental data since these alternative sampling and processing strategies open up new applications for magnetic resonance. In the future, these advanced techniques may also expand the use of low-field MR technology in process applications.

Relaxation time analysis
Spin-lattice (T 1 ) and spin-spin (T 2 ) relaxation time measurements have long been utilized to characterize porous materials (Brownstein andTarr 1977, Watson andChang 1997) and provide an in situ determination of pore-size distributions (Davies and Packer 1990). For several years, there has been extensive debate about the relative merits of selecting either T 1 or T 2 as a characterization of porous media (Watson and Chang 1997, Allen et al 1998, Aksnes and Gjerdåker 1999, Mitchell et al 2010a. In summary, it is easiest to differentiate between them on the basis of the fact that T 1 and T 2 processes are averaged over different temporal scales. T 1 probes variations occurring over milliseconds to seconds, whereas T 2 probes variations occurring on a microsecond timescale. For this reason, T 1 is often considered to be an indicator of the overall pore structure, while T 2 reflects local pore geometry and surface chemistry. Recently, two-dimensional (2D) relaxation time correlations have become popular for exploring surface interactions and transport at the microscopic scale in liquid-saturated porous materials. In particular, T 1 -T 2 and T 2 -T 2 correlations are finding widespread applications. In general, T 1 -T 2 correlations are used to characterize sorbate-surface interactions, while T 2 -T 2 correlations are useful for characterizing transport and exchange processes at the pore scale of microporous and mesoporous materials.
The T 1 -T 2 correlation  is of particular interest in chemical engineering because the ratio T 1 /T 2 yields information on the strength of surface interactions , Mitchell et al 2009a. In general, the stronger the interaction 4 (adsorption) between a liquid molecule and the pore surface, the larger the observed (averaged) T 1 /T 2 ratio of the liquid becomes. There are two models describing this trend: 1. The electron-proton interaction between surface-diffusing liquid molecules and paramagnetic species on the pore surface gives rise to different relaxation rates for the longitudinal and transverse processes. The ratio of the surface diffusion correlation time τ m and the surface residence time τ s is related to the T 1 /T 2 ratio, and therefore T 1 /T 2 becomes a measure of surface adsorption strength (Kleinberg et al 1994, Godefroy et al 2001. 2. In the absence of paramagnetic species on the pore surface, the ratio T 1 /T 2 is related to the degree of mobility of a surface-adsorbed molecule. T 1 and T 2 processes are sensitive to different frequencies of translational and rotational motion. Therefore, when the molecular motion is restricted by surface adsorption, the T 1 /T 2 ratio is altered (Liu et al 1991).
Overall, the exact value of the T 1 /T 2 ratio will be determined by a number of factors, including surface molecular mobility, concentration of pore surface paramagnetic species, partial saturation of pores, local pore geometry and bulk liquid properties. However, T 1 /T 2 does provide a qualitative comparison between similar samples (e.g. the same porous solid, different imbibed liquids), in which case it can be assumed that the reduction in T 1 and T 2 relaxation times due to local pore geometry and bulk liquid properties will be consistent and therefore will have a negligible influence on the ratio. However, it is important to ensure that the pores are saturated fully in these experiments (Allen et al 1998(Allen et al , 2001. In examples discussed here, the dominant parameters determining the T 1 /T 2 ratio are assumed to be surface molecular mobility and concentration of pore surface paramagnetic species, both of which are related to surface interaction strength, as discussed above. The T 1 -T 2 correlation measurement comprises a radio-frequency (rf) excitation pulse and T 1 recovery delay (Vold et al 1968) followed by a Carr-Purcell Meiboom-Gill (CPMG) echo train (Carr andPurcell 1954, Meiboom andGill 1958). The intensities of all the echoes in the CPMG portion of the sequence are recorded in a single shot. The T 1 recovery delay is incremented in successive experiments to construct a 2D data array. The overall experimental time is similar to that for a 1D T 1 measurement. It is usual to apply a 180 • excitation pulse and hence monitor the inversion recovery of the spin ensemble to determine the T 1 relaxation time dimension . However, this can be time consuming since inversion recovery requires the repeat delay t R = 5 × T 1 . For example, a typical T 1 -T 2 acquisition with 32 recovery times spaced logarithmically, a recycle delay of t R = 10 s and eight repeat scans will have a total experiment time of t exp ≈ 60 min. Saturation recovery can be utilized instead, whereby the excitation pulse is replaced by a comb of 90 • rf pulses. The length of t R then becomes unimportant. This allows a full T 1 -T 2 correlation with 16 recovery times spaced logarithmically, a recycle delay of t R = 1 s and two repeat scans to be acquired in t exp ≈ 3 min (Seland et al 2007), and allows dynamic systems to be studied (e.g. Song et al 2010). A requirement of the saturation recovery sequence is that the spin ensemble is well saturated in the x-y plane at the start of each acquisition. This can be achieved by ensuring that the 90 • pulse train contains at least 32 pulses separated by a variable time delay and variable amplitude homospoil gradients to prevent the formation of unwanted spin or gradient echoes. Once the 2D data set has been acquired, it is inverted using an appropriate numerical method. When the two data dimensions are separable (i.e. they do not share a common time base), the data can be inverted using the algorithm given by Venkataramanan et al (2002). The ratio T 1 /T 2 is then estimated from the correlated 2D distribution of T 1 and T 2 relaxation times. The spin ensemble is driven to some chosen equilibrium state during the DE portion of the sequence determined by the ratio τ 1 /τ 2 where τ 2 = 4τ SE and then the intensity of the spin echoes are recorded during the CPMG portion of the sequence.
An alternative, direct measurement of the T 1 /T 2 ratio is the driven equilibrium CPMG (DECPMG) sequence (Mitchell et al 2009a) shown in figure 1. An equilibrium magnetization is achieved by p repetitions of the driven equilibrium (DE) encoding portion of the sequence. The magnetization M DE following the DE encoding is where M 0 is the total (initial) magnetization and τ 1,2 are timing delays, as defined in figure 1. The equilibrium value can be checked by recording the echo generated during this portion of the sequence, although this is not required for the data analysis. The intensities of the CPMG spin echoes are recorded during the second portion of the pulse sequence. The amplitude of the qth echo is given by where t E = 2τ CPMG is the temporal separation of adjacent echoes (the so-called echo spacing). By comparison to the amplitude of a reference CPMG measurement where the magnetization is given by M C = M 0 exp{−qt E /T 2 }, the T 1 /T 2 ratio is determined as By inverting the DECPMG and reference CPMG measurements into T 2 distributions, equation (3) can be used to determine the T 1 /T 2 ratio as a function of T 2 . A DECPMG acquisition with eight repeat scans will have a total experiment time of the order of t exp ≈ 80 s, depending on the maximum T 1 and T 2 relaxation times being observed. As with saturation recovery, the recycle delay t R is unimportant because a sufficient number of loops p of the DE portion of the sequence will ensure that a predetermined equilibrium state is achieved regardless of the initial condition of the spin ensemble. However, this requires p ∼ 1000 typically and so the DECPMG experiment has a high rf duty cycle and so may not be suitable for all spectrometer hardware. Therefore, the T 1 -T 2 correlation with saturation recovery is still the T 2 -T 1 -δ pulse sequence. The T 2 dimension is encoded in a train of m spin echoes. This is followed by a double-shot T 1 sequence where n truncated FIDs are acquired after each small tip angle α rf pulse. The grey rectangle represents a homospoil gradient pulse. preferred generic pulse sequence for a rapid determination of the T 1 /T 2 ratio, particularly in high-field measurements; the DECPMG sequence offers an alternative in low-field applications where the rf power requirement per pulse is lower.
As mentioned above, a standard T 1 -T 2 acquisition has a similar experimental time to a 1D T 1 measurement. However, if chemical sensitivity is required, a full 3D T 1 -T 2 -δ will have a much longer experimental time because the T 1 and T 2 dimensions have to be incremented independently. This limitation is overcome in the T 2 -T 1 -δ experiment (Chandrasekera et al 2008). The pulse sequence is shown in figure 2 and comprises a CPMG echo train followed by a double-shot T 1 acquisition. The double-shot T 1 measurement provides a signal that decays exponentially with T 1 . A truncated free induction decay (FID) is acquired after a small tip angle α rf pulse. A homospoil gradient is applied to remove any remaining spin coherence prior to the next acquisition. This is akin to the Difftrain sequence described in section 2.2 (albeit without diffusion encoding gradients). If the FID is short, i.e. the sample has a broad spectral line as observed frequently in porous media, the acquisition window is sufficient to provide a chemical spectrum. The experimental time of the T 2 -T 1 -δ sequence is similar to a standard T 1 -T 2 correlation measurement with inversion recovery. For example, a T 2 -T 1 -δ acquisition with 32 CPMG loops, 32 T 1 recovery times, a recycle delay of t R = 10 s and eight repeat scans will have a total experimental time of t exp ≈ 70 min. The signal-to-noise ratio (SNR) of the T 2 -T 1 -δ data is obviously reduced due to the α pulses, although this sequence would be implemented typically on an intermediate-to high-field spectrometer system where the SNR is expected to be high.
Transport over microscopic length scales can be monitored using 2D relaxation experiments when liquid molecules move from one relaxation environment to another, such as between pores of different size. This motion is observed in T 1 -T 2 correlations in the form of 'exchange peaks' which exhibit a distinctly different T 1 /T 2 ratio compared to the rest of the system. These exchange peaks can be visualized more robustly in T 2 -T 2 exchange experiments (McDonald et al 2005, Washburn andCallaghan 2006). The T 2 -T 2 exchange experiment comprises two CPMG echo trains separated by an exchange (storage) time t S . During the exchange time, the spins are stored on the ±z-axis and restored as a stimulated echo. Therefore, 7 the spins undergo T 1 relaxation during the storage interval. A typical maximum exchange time in 1 H measurements is t S = 2 s, although this will depend on the T 1 of the sample. The number of echoes in the first CPMG train is incremented in successive acquisitions in order to construct a 2D data set; the amplitude of every echo in the second CPMG train is recorded. Following a 2D numerical inversion to form a correlated distribution of T A 2 against T B 2 , the majority of the signal lies on the diagonal where T A 2 = T B 2 . Exchange peaks are observed either side of this line. Utilizing a model for diffusive exchange in the presence of relaxation allows the exchange rate to be determined from the evolution of the exchange peaks as τ ex is increased in successive experiments (Monteilhet et al 2006). This model has been validated in water-saturated layered systems of packed beads where alternating layers were associated with different T 2 relaxation times due to differing surface chemistry of the beads (Mitchell et al 2007). The observed exchange rates corresponded to diffusion distances on the same scale as the width of the layers. A typical T 2 -T 2 acquisition, with 32 CPMG loops in the first dimension, a recycle delay of t R = 10 s and eight repeat scans, will have a total experiment time of t exp ≈ 60 min. The acquisition of a full 3D T 2 -T 2 -t S data set can be time consuming, although it is usual to acquire only a small number of different t S times. So, for example, a T 2 -T 2 -t S acquisition with eight exchange times will have a total experiment time of t exp ≈ 8 h.
One significant limitation of intermediate-and high-field relaxation time analysis of liquidsaturated porous media is the influence of internal gradients on the T 2 dimension (Mitchell et al 2010a). However, the internal gradients have no effect on T 1 relaxation time measurements (Washburn et al 2008). Internal gradients arise due to the magnetic susceptibility contrast between the liquid and solid, and scale with increasing field strength but inversely with pore size. As spins diffuse through the internal gradients, dephasing of the coherent spin ensemble occurs and this contributes to the loss in magnetization recorded during a CPMG echo train. The relaxation time estimated from such a decay is referred to as the effective transverse relaxation time T 2,eff . If the spins dephase by more than 2π radians by diffusing a length l g between successive echoes, then the resultant signal will be insensitive to details of the pore space and will be governed predominantly by diffusion and internal gradient strength. Depending on the pore size and internal gradient strength, there are three diffusion regimes that need to be considered: 1. The motional averaging (MAV) regime (Hürlimann et al 1995) is observed when the pore size l s is less than the dephasing path length l g and the diffusion path length l E associated with the echo time t E . The magnetization decay will be governed by an exponent that varies with l 4 s and t E . 2. The short time (ST) regime (Le Doussal and Sen 1992) is observed when the diffusion path length l E is less than the dephasing length scale l g and the pore size l s . The magnetization decay will be governed by a surface relaxation T 2 exponent that varies with t E and a diffusive exponent that varies with t 3 E . This is the optimum regime for T 2 measurements because the diffusive and surface relaxation exponents vary with different powers of t E and can therefore be distinguished (Mitchell et al 2010b). 3. The pre-asymptotic localization (LOC) regime (de Swiet and Sen 1994) is observed in the presence of large internal gradients where the dephasing path length l g is shorter than either the diffusion path length l E or the pore size l s . The magnetization decay is governed by a diffusion exponent that varies with t x E where 1 < x < 3 and a surface relaxation exponent term that varies with t E . As such, it is still possible, under certain conditions, to separate 8 the surface relaxation and diffusion components in the pre-asymptotic LOC regime. The asymptotic limit of the LOC regime is not observed generally in experimental studies of 2D and 3D pores due to the rapid loss of magnetization.
If the length scales l E , l s and l g can be estimated from the experimental parameters, it is possible to determine in which regime diffusion will occur. However, the boundaries between the regimes are not well defined because the sample is likely to contain a distribution of pore sizes l s and effective internal gradients g eff . We note that the presence of strong internal gradients can lead to the observation of multiple effective relaxation environments within individual pores. This is most pronounced in high-field T 2 -T 2 experiments where exchange between regions of different internal gradient strengths is observed. It has been shown that this exchange rate can be used as an alternative measure of pore size (Mitchell et al 2010a).
If diffusion is occurring in or near to the ST regime, then it is possible to determine the internal gradient distribution and hence correct any T 2 distribution (e.g. as found in the second dimension of a T 1 -T 2 correlation) for the effects of internal gradients. The required data are simply a set of CPMG decays, each with a different echo time t E . The acquisition of such data with 32 different echo times, a recycle delay of t R = 10 s and 8 repeat scans will have an experimental time of t exp ≈ 60 min. The amplitude of the CPMG decays is recorded at a fixed time t F = nt E , where both n and t E vary to ensure that t F is constant. The surface relaxation T 2 exponent varies with t E , so the loss of magnetization due to T 2 relaxation is constant at t F , regardless of t E . However, the diffusion exponent varies with t 3 E , so the loss of magnetization due to diffusion at t F will vary with t E . Therefore, the signal intensity at t F (obtained from the 32 CPMG decays) can be fitted to generate an internal gradient distribution (Sun and Dunn 2002). The magnetization decay due to diffusion can then be determined for an arbitrary t E and deconvolved from a CPMG decay to leave only the T 2 contribution from surface relaxation. This is subsequently inverted to provide the 'true' T 2 relaxation time distribution. Alternatively, an appropriate 2D numerical transform can be used to generate directly a correlation between internal gradient strength g eff and 'true' T 2 distributions (Mitchell et al 2010b). This correction on measurements of T 2 distributions allows porous media samples to be studied at higher magnetic field strengths than considered possible previously, and is important for determining the correct T 1 /T 2 ratio in intermediate-to high-field experiments. Examples of applying T 1 -T 2 correlations to obtain robust T 1 /T 2 ratios in applications to materials processing and catalysis are discussed in section 3.1.

Flow propagators
The NMR flow propagator is a probability distribution of molecular displacements measured as a function of time (Kärger and Heink 1983). Flow propagators are measured using conventional pulsed field gradient (PFG) diffusion measurements; typically, the alternating pulse gradient stimulated echo (APGSTE) sequence (Cotts et al 1989) is used, as shown in figure 3. A pair of bipolar gradients encodes the spin phase during the first spin echo portion of the sequence. The spins are then stored on the ±z axis for a time T and restored to form a stimulated echo. The phase is decoded by another pair of bipolar gradients. The advective displacement that occurs over the observation time is encoded in the phase shift of the spins and detected as a reduction in the absolute signal intensity. By incrementing the encoding gradients, the acquisition steps through q-space. The flow propagator is obtained from the Fourier transform of the q-space data. The bipolar gradients are preferential in porous media studies as they reduce the influence of background magnetic field inhomogeneities.
There are three main limitations of the APGSTE sequence as applied to studies of porous media: (i) the maximum observation time, (ii) the experimental duration and (iii) quantitative signal acquisition. Different solutions are available to overcome each of these limitations. These will now be considered in turn.

Maximum observation time.
In general, we wish to study the evolution of molecular transport within a porous structure for as long a time as possible. For flow propagator measurements of 1 H nuclei, the maximum observation time max is limited to approximately 4 × T 1 . Even in bulk liquids, 1 H T 1 times are rarely longer than 2-3 s, so max ≈ 8-12 s. In porous media, T 1 times are often reduced, limiting max further. This can be overcome by measuring 13 C flow propagators: T 1 times of 13 C nuclei can be significantly longer than those of 1 H, allowing the maximum observation time to be extended up to max ≈ 30-35 s. However, the 1.07% natural abundance of 13 C, along with the reduced Larmor frequency, means that the MR sensitivity is 5870 times less than for 1 H. In many applications, it is impractical, physically or economically, to use 13 C isotopically enriched samples. In such cases, distortionless enhancement by polarization transfer (DEPT) (Doddrell et al 1982) can be employed to boost the signal from the natural abundance of 13 C. DEPT has been incorporated into the flow propagator measurement to facilitate detection of the carbon signal (Akpa et al 2007). This requires the addition of a polarization transfer pulse sequence prior to the APGSTE 13 C experiment, as shown in figure 4.

Experimental duration.
APGSTE experiments can be time consuming, especially if long observation times and a large number of q-space points are to be explored. For example, the acquisition of eight flow propagators each containing 128 q-space points, with a recycle delay of t R = 3 s, a maximum observation time of max = 2 s and eight repeat scans, requires t exp ≈ 8 h. One option for acquiring multiple flow propagators rapidly is the Difftrain pulse sequence (Stamps et al 2001). The Difftrain sequence allows a single q-space point to be acquired at multiple observation times following a single excitation. Although this does not speed up the acquisition of one propagator, it does allow 8 or 16 propagators to be acquired in the same experimental time as a standard APGSTE measurement. So the Difftrain acquisition of the eight flow propagators described in the previous example has a total experimental time of only t exp ≈ 70 min. The various implementations of Difftrain have been reviewed recently by Mitchell and Johns (2009); in figure 5, we show only the most recent and universally applicable version. Superficially, the Difftrain sequence looks similar to the APGSTE sequence, except that the stimulated echo restore rf pulse now provides a small tip angle of, say, α ∼ 10 • . The decoding portion of the APGSTE sequence is repeated n times and the observation time of each acquisition sums as where T start and T inc,n are defined in figure 5. The inclusion of additional homospoil gradients around the 180 • rf pulses improves the phase coherence of the detected signal (Buckley et al 2003). During each n repetition, the amplitude of the α rf pulse and homospoil gradients can be varied to optimize the acquisition parameters. The total observation time is given by equation (4).
It is possible to adjust the amplitude of the α rf pulse in each of the loops. If the T 1 of the sample is known, the values of α n can be determined from to provide a uniform SNR across all the propagators by increasing the amplitude of the signal acquired at longer observation times (Davies et al 2007).
The speed of the propagator measurement can also be increased greatly by under-sampling q-space. If the full probability distribution is not required, the first, second and third order moments (approximating to the mean ζ , rms width σ and adjusted skewness γ of the displacement distribution) can be obtained from a small selection of data with low q-space values (Scheven et al 2005). Using this approach, the moments of eight flow propagators can be determined from 16 q-space points acquired using the Difftrain pulse sequence in t exp ≈ 9 min (Mitchell et al 2008a). The moments are obtained from the phase ϕ and magnitude of the signal |S| as and ln |S(q)| = − 1 2 σ 2 q 2 .
Therefore it is only necessary to sample a range of q-space that is well fitted by the cumulative expansion of q n up to n = 3 (Mitra et al 1992).

Quantitation.
The propagator measurement is inherently non-quantitative due to relaxation weighting. T 2 relaxation occurs during the spin echo portions of the pulse sequence, when the phase encode/decode gradients are applied, and T 1 relaxation occurs during the stimulated echo storage time T, as defined in figure 3. Removing the relaxation weighting is not straightforward since the relaxation times vary with flow rate in porous media: stagnant (i.e. stationary) spins will undergo more pore-surface interactions than flowing spins, and conversely flowing spins will experience more susceptibility-induced magnetic field gradients. An analytic correction has been described by Scheven et al (2005). Firstly, the ratio θ = ζ / ζ 0 is determined where ζ 0 is the expected mean displacement estimated from the pore velocity. The entire probability distribution is rescaled by a factor θ −1 and a Dirac delta function δ dd of amplitude (1 − θ −1 ) is added to account for lost signal from stagnant spins. The corrected probability distribution P (ζ, ) is given by and the corrected second σ and third γ moments can be determined accordingly.
More recently, an empirical correction has been demonstrated in which T 1 and T 2 are determined as a function of displacement and each data point in the propagator is then corrected independently (Mitchell et al 2008b). This is an improvement over the analytic correction described above because it makes no assumption about the uniformity of the correction across the probability distribution. The determination of a ζ -T 2 correlation is straightforward and requires only the addition of a CPMG echo train at the end of the APGSTE acquisition (Britton et al 2001(Britton et al , 2004. However, it is advisable to include an additional stimulated echo prior to the CPMG sequence to remove unwanted spin coherences (see figure 6). The time base of the CPMG acquisition is adjusted to account for all transverse relaxation that occurs during the pulse sequence. Therefore, the intensity of the first echo is acquired at an effective time 4τ + 2τ CPMG , as defined in figure 6. By applying a suitable fitting algorithm to the CPMG decay, the signal intensity without T 2 weighting is obtained. This is equivalent to the projection onto the ζ axis of a ζ -T 2 correlation determined with the modified CPMG time base.
Obtaining a flow propagator without T 1 weighting is less straightforward. Due to the stimulated echo inherent in the APGSTE sequence, T 1 recovery times in a T 1 -ζ experiment are limited by the stimulated echo duration T ; notwithstanding, an appropriate pulse sequence has been presented in the literature (Mitchell et al 2008a). Instead of measuring a T 1 -ζ correlation directly, the ratio T 1 /T 2 for the sample is determined as a function of T 2 using an appropriate method such as a T 1 -T 2 correlation or a DECPMG experiment, as described in section 2.1. Once the T 1 /T 2 ratio has been determined, it is used to adjust the amplitude in a ζ -T 2 correlation to account for longitudinal relaxation during the time T + d in the ζ -T 2 pulse sequence (see figure 6). The projection onto the ζ axis of the rescaled ζ -T 2 correlation will then provide a probability distribution that does not contain T 1 or T 2 relaxation weighting. A similar process could be used to correct just the low q-space points in order to obtain corrected moments using the method described above. However, unlike the analytic correction, this empirical technique is not compatible with the rapid Difftrain acquisition.

Imaging
The theory and methods of MRI are well documented in various texts (see e.g. Callaghan 1991, Haacke et al 1999. The majority of advances in MRI have occurred in medical imaging applications. These pulse sequences can be used directly in studies of porous media, although the different requirements of chemical engineering have ensured that significant modifications have been developed to gain more information from these imaging sequences. For example, medical images often lack quantitative signal intensities, relying instead on signal contrast due to relaxation or diffusion phenomena for interpretation. This is considered insufficient generally in studies of mass transport where quantification is of primary interest. In a typical MRI experiment, the voxel size is limited to 50-100 µm approximately in any dimension, and therefore direct imaging of structure in microporous or mesoporous networks is not possible. Specialized imaging equipment (magnets and gradient coils) has been developed to provide spatial resolution down to the order of a few microns in uniform samples, as reviewed by Mitchell et al (2006). However, the magnetic susceptibility contrast inherent in porous materials is often the limiting factor in high-resolution imaging experiments. In general, studies of microporous and mesoporous systems rely on direct spin density images of the imbibed fluid within the pore space to map macroscopic variations in voidage. By overlaying equivalent, spatially resolved maps of T 1 , T 2 or molecular self-diffusion coefficient D, information about structure-transport relationships within the system can be deduced. Clearly, when considering macroporous systems, the spatial resolution available makes direct imaging of the pore network structure possible and examples of this will be discussed in section 3.2.
2.3.1. Imaging structure and texture. Images of the internal structure of a porous sample, by which we mean spin density maps of a liquid (or gas) imbibed within the pores, can be obtained using the basic spin-warp pulse sequence for a full 3D image acquisition, as shown in figure 7 (Edelstein et al 1980, Johnson et al 1983. The spins are excited initially by a 90 • rf pulse, after which imaging gradients of amplitude G read and G phase are applied in the read and twophase encode directions, respectively. Once the spins have evolved for half the echo time t E /2, an 180 • rf refocusing pulse is applied to form a spin echo. A second read gradient is applied during the acquisition of the spin echo, centred at time t E after the rf excitation pulse. The signal intensity across the spin echo is modulated by the spectral line width (1/T * 2 ) and the total signal intensity is weighted by T 2 . Ideally, the recycle delay t R = 5 × T 1 if the maximum signal (magnetization) is to be obtained. However, a reduced t R may be used to provide a compromise between maximum signal and minimum acquisition time (e.g. t R = 3 × T 1 ); in this case, the echo amplitude will be T 1 weighted as well.
The Larmor frequency ω of the spins at position r is related to the applied magnetic field gradient vector G by where γ is the gyromagnetic ratio of the resonant spins and B 0 is the magnitude of the static magnetic field strength. The frequency domain spatial image ρ(r) is obtained by the 3D Fourier transform of the time domain k-space data S(k) such that The reciprocal space vector k is defined as where t p is the gradient duration (encoding time). The spin echo, acquired in the presence of the read gradient, is a single line through one dimension of k-space. A large number of data points are recorded more readily in the read direction due to the nature of the spin echo acquisition. Consequently, this gradient direction is usually chosen to correspond to the longest sample dimension. In order to sample the other two dimensions, it is necessary to increment the two-phase encode gradients independently and sequentially from −G phase to +G phase . Therefore, the acquisition of a full 3D data set in this manner can be time consuming, depending on the resolution required. For example, if the image were to contain 128 × 64 × 64 voxels with a recycle delay of t R = 10 s, and 1 scan, then the total experimental time would be t exp ≈ 11 h.
Although it can take a long time to acquire a full 3D image using the spin-warp experiment, relative (voxel-to-voxel) quantification is retained because the images have a uniform relaxation contrast.
If quantification (absolute signal intensity) is important, then pure phase encoding pulse sequences based on single point imaging (SPI) techniques can provide a solution (Emid and Creyghton 1985). Here, a single data point is acquired with a minimum encoding time t p , ensuring that the spin density in each image voxel is not weighted significantly by T * 2 , T 2 or T 1 relaxation times. However, a 3D SPI image of the same resolution as the spin-warp imaging example given above would have a total acquisition time of t ex ≈ 8 1/2 weeks. Rapid versions of the SPI pulse sequence have been devised, such as single point ramped imaging with T 1 enhancement (SPRITE) (Balcom et al 1996), which allow a 3D image to be acquired in a matter of minutes. The SPRITE class of sequences is used typically to acquire lower resolution images (e.g. 32 × 32 × 32 voxels) and relies on the sample having a short T 1 relaxation time of the order of T 1 ∼ 10-100 ms. Additionally, the gradient duty cycle of a 3D SPRITE acquisition is high, preventing implementation on many commercial imaging magnets.
If quantification is not a primary concern in a structural image, rapid imaging techniques can be employed as an alternative to the spin-warp or SPI-type sequences. One such method is the rapid acquisition with relaxation enhancement (RARE) pulse sequence (Hennig et al 1986), shown in figure 8 (also referred to as turbo spin echo imaging). Shaped rf pulses are applied during the slice gradients to select a volume slice of the sample. A line of k-space is then recorded by applying a read gradient during the acquisition of a spin echo. The position of this line is determined by the phase encoding gradient, G phase , applied prior to the echo. The phase encoding is subsequently removed by applying a phase decoding gradient of amplitude −G phase after the echo. A train of echoes is generated by the repeated application of shaped 180 • rf pulses, akin to the CPMG echo train. The amplitude G phase is incremented with each echo to step through r k-space lines during the acquisition of the echo train. In figure 8, four repetitions of the echo acquisition, denoted by loop A, are shown. This acquisition is therefore said to have a RARE factor of r = 4. As the echo train evolves, the signal amplitude will decrease inherently due to T 2 relaxation, and so the acquired image will contain T 2 contrast. Since the maximum signal is acquired at the centre of k-space, it may be advantageous to employ two echo trains where the phase encode gradients range from 0 → −G phase and 0 → +G phase , ensuring minimum T 2 weighting on the central k-space data point (equivalent to total image intensity). If the T 2 relaxation time is long, large RARE factors can be employed, allowing the rapid acquisition of high resolution images. For example, a 128 × 64 × 64 voxel image, employing a RARE factor of r = 64, an echo time of t E = 3 ms, a recycle delay of t R = 1 s and 1 scan, will have a total acquisition time of t exp = 76 s. Due to the use of a slice-selective excitation pulse, the acquisition of the different slices can be interleaved, so it is not necessary to wait for t R = 5 × T 1 between successive slice acquisitions. Therefore, the t R used in such a RARE acquisition is determined only by the maximum duty cycle of the gradient coils and amplifiers. However, in porous media applications, it may be necessary to reduce the RARE factor in order to accommodate the short effective transverse relaxation times T 2,eff encountered commonly in these samples.
A modified version of the RARE sequence allows the acquisition of multiple images following a single rf excitation pulse. This is the single excitation multiple image RARE For the SEMI-RARE sequence, A is repeated n × r times. Every r repetition of A, the G phase increments are repeated to acquire n 2D images.
(SEMI-RARE) pulse sequence , which is shown in figure 8. In the SEMI-RARE implementation, n images are obtained, each containing r pixels in the phase direction, in a single echo train. The sequence of phase gradient increments required for a single image-slice acquisition is repeated multiple times during the evolution of the spin ensemble. Like the RARE sequence, SEMI-RARE is subject to relaxation contrast, except that the T 2,eff weighting will be carried across multiple images. However, the increased speed of the acquisition allows rapid time-of-flight (ToF) imaging of dynamic systems that could not be resolved using the standard RARE sequence due to the recycle delay between successive acquisitions of the same image slice. SEMI-RARE was used to image two-phase flow in parallel channel ceramic monoliths (Mantle et al 2002). In later work, SEMI-RARE was extended such that liquid velocities up to 5-10 m s −1 in gas-liquid two-phase flow could be tracked within these monoliths .
As well as being optimized for ultra-fast image acquisition, the standard RARE sequence has also been adapted for other specialized applications, an example being the rotationally compensated RARE (ROTACOR) pulse sequence (Sederman et al 2004a). During the ROTACOR sequence, the gradient co-ordinates are rotated to allow undistorted images of rotating liquid samples to be acquired. This approach was first applied to image the deformation of water droplets in silicone oil during shear within a wide-gap Couette cell. Accurate imaging of the droplet shape as a function of shear rate enabled the in situ determination of the interfacial tension between the two fluids.
Another rapid imaging sequence used widely in medical applications is echo planar imaging (EPI) (Mansfield 1977). This comprises a single slice-selective rf excitation pulse followed by a train of alternating read gradients with blipped phase increment gradients between each read gradient. A 128 × 64 × 64 voxel image, with a recycle delay of t R = 1 s and 1 scan, will have a total acquisition time of less than t exp = 70 s. However, as a train of gradient echoes are being acquired, the final image is T * 2 weighted. The broad spectral lines encountered in  Figure 9. The PEPI-hybrid imaging sequence allows the acquisition of EPI images in short T * 2 materials. Signal acquisition begins after ACQ; the first (unobserved) spin echo is not shown. The acquisition sequence is repeated n/2 times with a blipped gradient applied on the G phase axis to step through k-space.
porous media therefore result in distortions and blurring in EPI images. The π -EPI (PEPI) sequence has been demonstrated to provide improved performance in application to short T * 2 materials (Guilfoyle et al 1992, Peters et al 1996. Conventional EPI utilizes gradient echoes to scan rapidly through k-space, whereas PEPI uses spin echoes generated by a series of 180 • rf pulses, combined with gradient reversal, to step through k-space. An additional modification is the PEPI-hybrid sequence in which the gradient reversal is removed (Guilfoyle et al 1996; see figure 9). This pulse sequence is very similar to RARE in its implementation and total acquisition time. It does, however, rely on good-quality hard 180 • rf pulses to prevent the accumulation of unwanted phase coherences leading to artefacts in the final image. The PEPIhybrid sequence has been used to acquire images of porous materials with a resolution better than 100 µm on a standard spectrometer (Manz et al 1999b).
We note here that RARE is much more robust (i.e. artefact free) in application to systems characterized by heterogeneous magnetic susceptibility typical of those found in chemical engineering. Furthermore, the decay of the MR signal in a RARE experiment is dependent on T 2 (and, of course, to some extent on T 1 ) and not T * 2 , as is the case for EPI. In the magnetically heterogeneous systems of interest to us, T * 2 is considerably shorter than T 2 , and therefore RARE allows acquisition of data of higher SNR over longer timescales. This is important if one wishes to investigate the time evolution of a system. The differences between these pulse sequences and hence their relative strengths are discussed in more detail elsewhere Sederman 2003, Gladden et al 2006).
EPI-based pulse sequences represent the fastest data acquisitions available if we are to sample the whole k-space array. However, it is well known that only a sub-set of the data points in k-space contains the majority of the signal intensity and therefore most of the information in the image obtained after Fourier transformation. CS allows accurate reconstruction of an under-sampled image with prior knowledge that the image is sparse when transformed to some other domain, e.g. by a wavelet or curvelet transform (Candès and Guo 2002). This is a reasonable assumption for k-space MRI acquisitions. If the image to be reconstructed is stacked as a vector x, and is the transform operator, then the reconstruction is obtained by solving the optimization problem min x 1 where the 1 -norm x 1 = i |x i | is a proxy for sparsity. This optimization is subject to the constraint Fx − y 2 < ε where F is the under-sampled Fourier transform operator, y is the acquired k-space data (as a vector) and ε is a threshold equal to the expected noise. The image can be determined using a variety of methods, including iterative reconstruction (Fornasier and Rauhut 2008) and projected conjugate gradients (Lustig et al 2007). A first example of CS being applied to real MR data reconstruction was demonstrated in conjunction with SPI for mapping moisture distribution in a porous food wafer (Parasoglou et al 2009). Each image was acquired in 13 min; an equivalent full k-space sampled image acquisition would have taken 39 min.
In studies of transport in porous media, it is often required that different chemical species be monitored simultaneously. This can be achieved by combining MR spectroscopy with MRI protocols in chemical shift imaging (CSI) (Maudsley et al 1983). Although spectral contrast can be achieved by including a frequency-selective shaped rf pulse in any standard imaging sequence, true CSI can only be achieved by acquiring a FID in the absence of any imaging gradients. In this way, a spectrum will be associated with each pixel or voxel in an image. The easiest implementation of CSI involves incremental phase encoding to step through k-space, as demonstrated in figure 10 for a 1D chemical shift profile. This can be extended up to three dimensions by the addition of two more orthogonal G phase gradients. However, this type of acquisition is time consuming (equivalent to a full 3D SPI image, as discussed above) and alternative methods are available to reduce the experimental time.
One alternative to three incremental phase encode gradients is the use of a slice-selective gradient and a shaped rf pulse to provide a series of 2D chemical shift image slices. This would allow a 3D CSI with 32 × 32 × 32 voxels, a FID acquisition time of t acq = 50 ms, a recycle delay of t R = 800 ms and 1 scan to be acquired in t exp ≈ 8 h. A further extension of this concept is the use of slice-selective gradients on all three axes to provide volume-selective spectroscopy  Figure 11. Volume-selective CSI pulse sequence. Three soft rf pulses are applied, all with corresponding slice-selective gradients to select a voxel of interest. A full FID is acquired to provide the chemical spectrum within that voxel. (Kimmich and Hoepfel 1987); an example of such a pulse sequence is shown in figure 11. This allows spectra to be acquired rapidly from selected regions of interest within a sample, such as when mapping chemical composition in a fixed-bed reactor during a reaction (Yuen et al 2002). Even when accumulating several repeat scans to increase the SNR, the spectrum from a single voxel can be acquired in approximately t exp ≈ 3 min.
The localized acquisition inherent in CSI can be advantageous in porous media studies since the local magnetic field variations will be less than the macroscopic magnetic field variations across the entire sample. Therefore, the spectrum from a single CSI voxel may be better resolved (even though the SNR is lower per scan) than a spectrum from the bulk sample since the individual lines will be narrower. Examples of localized acquisition overcoming spectral broadening due to magnetic field heterogeneity have been reported by Halse and Callaghan (2007) and Lim et al (2011).
It is important to note, however, that CSI techniques will not overcome line width limitations imposed by local (i.e. pore-size scale) magnetic field heterogeneities. Therefore, even with localized acquisition, a 1 H spectrum may be broadened significantly due to magnetic susceptibility contrast within individual pores, preventing the resolution of individual spectral features. Hence 13 C spectra offer a good alternative since the spectral lines of carbon exhibit a much larger chemical shift range and individual lines are therefore resolved easily even in the presence of susceptibility-induced internal gradients. Spatially resolved 13 C DEPT can be achieved by acquiring a 13 C FID in the presence of a broadband 1 H-decoupling rf pulse following the DEPT sequence and phase encode gradients (see figure 12; Yeung and Swanson 1989). In chemical mapping in catalytic reactors, for example, the use of polarization enhancement is essential because the large volumes of feed material make isotopic enrichment prohibitively expensive.

Velocity imaging.
Macroscopic transport within a porous structure can be visualized using velocity imaging. Velocity mapping is based on the concept of dynamic MR microscopy (Callaghan and Xia 1991).  Figure 13. Velocity-encoded MRI pulse sequence. Each pixel in a slice is attributed a characteristic velocity by incrementing the velocity encoding gradients g * . These are applied on the gradient axis corresponding to the direction of flow. The gradient timings δ and , and maximum gradient amplitude g max , are chosen to span a suitable range of q-space. The first echo occurs at time t E = 2τ and the second echo is acquired at time t E = 2τ following a slice-selective soft 180 • rf pulse.
The velocity imaging sequence in figure 13 (Sederman et al 1997) provides a 2D image slice generated by the G read , G phase and G slice gradients. One component of the velocity vector v is encoded and decoded by the inclusion of additional gradients g * placed on either side of a 180 • rf refocusing pulse. The g * gradients have to be applied separately on all the three axes to acquire the full velocity vector. These gradient pulses are of duration δ and separated by the observation time . The experiment proceeds by incrementing the g * steps through the q-space domain (Callaghan 1991). The Fourier transform of the resulting data yields the flow velocity distribution. The peak in the Fourier transformed data occurs at where N is the number of digital data points in the Fourier transformed velocity dimension, g max is the maximum value of g * and n D is the number of g * increments. From equation (12), the modal velocity component v of each pixel can be determined. 3D velocity maps can be obtained either by acquiring multiple slices or by replacing G slice with a second incremental G phase channel. The velocity v of fluid flowing in each voxel is obtained by comparing two images acquired with different encoding gradients g 1 and g 2 . The phase shift P s between the two images, and hence a measure of average velocity, is given as (Sederman and Gladden 2001a) This method allows the acquisition of 3D velocity images in the minimum experimental time. The total acquisition time will be determined by the number of pixels/voxels required, the recycle delay (as in standard imaging sequences), and also the observation time . For example, a 3D flow map containing 128 × 64 × 64 voxels, with a recycle delay of t R = 1 s, = 15 ms and two velocity encode steps, will have a total acquisition time of t exp ≈ 3 h. However, such a velocity map will still contain only one component of the velocity vector v. To obtain the velocity vector v = v x + v y + v z , three separate velocity maps need to be acquired, so t exp ≈ 9 h.
An extension to the basic velocity imaging is the gradient echo rapid velocity and acceleration imaging sequence (GERVAIS) (Sederman et al 2004b). The pulse sequence is shown in figure 14 and comprises an imaging section based on Blipped EPI, also known as MBEST-EPI (Howseman et al 1988), preceded by the velocity encoding. This portion of the sequence is repeated after the initial excitation to allow the determination of spatially resolved velocity autocorrelation functions. As such, this allows the acquisition of all three velocity components in an image slice following a single rf excitation.
The GERVAIS experiment is ideally suited to probing velocity fluctuations in turbulent systems with low Reynolds numbers, e.g. Re < 5000. Such systems include single-phase flow through a fixed-bed reactor (Sains et al 2005). In the original GERVAIS implementation, it is necessary to acquire a zero velocity (no flow) reference data set due to phase shifts inherent in the pulse sequence. However, this limitation has been overcome in the snap-shot GERVAIS (ssG) pulse sequence that acquires all five images required to reconstruct the three velocity components following a single excitation (Tayler et al 2010). The ssG sequence allows a velocity map to be acquired in a total experimental time of t exp ≈ 125 ms. This pulse sequence has been applied to acquire snapshot velocity fields within a moving oil droplet.
A further recent development has been the implementation of a snapshot velocity imaging pulse sequence based on the RARE pulse sequence. This approach, although not as fast as an EPI-based measurement, has all the advantages of robustness in implementation for chemical engineering applications of RARE relative to EPI that were discussed earlier. This fast RARE velocity pulse sequence has been demonstrated in application to imaging in situ the flow The GERVAIS pulse sequence. A 2D slice-selective image is generated by the G read , G phase and G slice gradients, based on a blipped EPI pulse sequence. Velocity gradients g * are applied in pairs around 180 • rf refocusing pulses on whichever gradient axis corresponds to the direction of motion. The EPI/velocity portion of the sequence is repeated n − 1 times following the initial excitation. This repeating unit has duration τ m and a set of gradient echoes are acquired in the time interval τ i .
field around dissolving controlled-release pharmaceutical delivery systems in a commercial dissolution test environment (Shiko et al 2010). Finally, we note that CS has also been implemented successfully for the first time for velocity imaging. The principle of the method is the same as described earlier for applying CS in normal spin density image acquisitions. However, further development was required to implement the method for the phase data acquired in velocity imaging. Sampling of only 30% of the k-space raster produced images which are indistinguishable from those calculated from fully sampled data arrays (Holland et al 2010a). The concept of under-sampling is of particular interest to us currently and within this activity we see significant potential in exploring the use of Bayesian approaches in MR. A first application of Bayesian methodology has been applied to characterizing bubble-size distributions in two-phase gas-liquid bubbly flows. The potential decrease in data acquisition times is truly significant here. Instead of acquiring a 256 × 256 data array and analysing the subsequent image, the bubble-size distribution is determined from acquisition of just 40 k-space data points (Holland et al 2010b). The Bayesian approach provides a sampling strategy that allows us to measure only the data we need to test a hypothesis and as such means that a full image acquisition is no longer required. In the current application, we simply ask the question: are our chosen k-space data points consistent with a proposed bubble-size distribution? This provides a further advantage because the system is no longer required to be in steady state over the timescale of an image acquisition; only the bubble-size distribution itself needs to remain constant during the measurement. For the system we have studied most recently, this new approach enabled us to determine bubble-size distributions in gas-liquid bubbly flows up to gas fractions in excess of 15%, significantly higher than can be studied by optical or EPIbased MR techniques.

Applications
For the purposes of this review, we separate the applications into two broad categories: (i) microporous and mesoporous systems and (ii) macroporous systems. In the former category, we consider recent insights into improving the MR characterization of transport in rock cores. We will also describe how MR studies of transport in model porous media, such as bead packs of small particles (<1 mm typically), can be used in combination with numerical simulations to enable us to not only validate the simulation codes themselves but also to use the codes to probe the dominant influences on transport in porous media over a hierarchy of length scales. Bead packs have often been chosen as a model porous medium of relevance to rock cores since they enable us to study a chemically 'clean' system where the inter-relationships of microstructure and transport are more readily investigated. As soon as real materials are considered, spatial variation in impurities, in addition to the issues associated with short relaxation times resulting from paramagnetic species and internal magnetic field gradients, make the studies of the fundamentals of structure-transport relationships far more challenging. We will also briefly discuss how the MR techniques developed to study rocks are now being used more widely in the optimization of manufacturing processes for building materials, and in understanding catalyst behaviour. In the overview of macroporous materials, we will consider porous media comprising packings of large elements (>1 mm typically), which are used to represent researchscale fixed-bed catalytic reactors. We will include a brief synopsis of work on fluidized granular media. The first example we present is from the third type of study, which also builds on the earlier discussion of dispersion in porous media. Manz et al (1999a) investigated the transport of water within unconsolidated packings of glass spheres such that a range of Peclet numbers (Pe = 40-3310) could be studied. The dispersion behaviour was then correlated with the radial distribution functions of the pore space calculated directly from 1 H spin density images of the 1 mm diameter beads in the pack. Dispersion occurring transverse to the direction of superficial flow was found to be dominated by Taylor dispersion. In the axial (flow) direction, non-Fickian effects were seen to become important, even for flow within a single pore. The axial dispersion coefficient approached an asymptotic limit characterized by a length scale that corresponded to the pore-to-pore packing correlation length in the direction of flow, suggesting that the non-Fickian dispersion reached a limiting behaviour after just one mixing between streams merging from two interconnected pores. This work was extended by comparing the predictions of a lattice Boltzmann code with MR velocity images of the flow field within a bead pack of 1 mm diameter spheres. Lattice Boltzmann was selected for this study because we were considering slow flows in a complex pore geometry-the fact that lattice Boltzmann is implemented without the need to create a mesh for the pore space is a great benefit in such applications. By importing the MR image of the porous medium directly into the flow simulator, a direct comparison of the predicted and measured flow fields could be made; the lattice Boltzmann code was also extended to predict the MR propagators associated with its simulated flow field. Given the excellent agreement between experimental data and the predictions of the lattice Boltzmann code, it was possible to use different implementations of the code to explore the origin of the dispersion over different length scales within the pore structure. By introducing normalized parameters, hydrodynamic dispersion over normalized displacement length scales spanning nearly 4 orders of magnitude (182 < Pe < 350; 0.4 < Re < 0.77) was investigated. Mechanical dispersion was shown to dominate over most of the range with holdup playing a significant role in the hydrodynamics at the largest length scales investigated.

Microporous and mesoporous
Two illustrations of how MR can be used to provide input parameters to simulations are given in the following examples: 1. A combination of MR and lattice Boltzmann was used to understand how a polymer treatment fluid deposits within a rock core structure, thereby reducing the permeability of the rock. This mechanism is of relevance to optimizing the treatment fluids formulated for use in enhanced (chemical) oil recovery (EOR) processes. To optimize the EOR, it is necessary to determine where the polymer deposits in the pore structure and how it consequently modifies the permeability of the rock. Two extreme scenarios were considered: (i) uniform deposition of the polymer on the pore surface and (ii) plugging of pores associated with stagnant flow. By comparing the MR propagators predicted according to each of the two scenarios with the observed MR propagators from the treated rock, it was possible to identify the most likely deposition mechanism-in this case, the experimental data were consistent with pore-plugging . 2. A different modelling strategy was employed to utilize 3D MRI images in both establishing a valid simulation lattice and identifying the essential details of the implementation of a cellular automata code. This code was used to describe the dissolution of hydrocarbon ganglia trapped within a water-saturated porous medium (Johns and Gladden 2002).
Even when the spatial resolution of an image is such that individual pores cannot be resolved, useful information can still be obtained from spin density images and T 1 maps of porous media; the T 1 being related to the average pore size (Brownstein and Tarr 1977) within a given pixel. A typical porous medium which presents this type of challenge is a porous oxide catalyst. Gladden et al (1995) showed that the T 1 maps of given batches of catalyst were characterized by a specific fractal dimension. Random walks performed on simulation lattices, constructed with the same fractal dimensions that characterize the catalysts, predicted molecular diffusion behaviour consistent with PFG NMR measurements of transport within the same materials . A similar approach, using the spin density image directly as a simulation lattice, has also been explored (Rigby and Gladden 1998). Such representations of meso-scale variation in porosity provide simulation lattices from which predicted transport properties are in good agreement with MR measurements of tortuosity and transport. This observation suggests that we do not necessarily have to be able to characterize a porous medium at the detailed level of the structure of individual pores to predict its in-use behaviour.

Transport in rocks.
This section addresses a sub-set of ongoing research in the area of microporous materials. Here we will consider primarily examples of obtaining quantitative information about the microstructure and transport characteristics of rock cores using relaxometry and propagator measurements. First, let us consider the acquisition of T 2 data and how these data are influenced by magnetic field strength. In section 2.1, it was noted that T 2 was influenced by internal magnetic field gradients that scale with magnetic field strength. Therefore, we ask the obvious question: why do we use T 2 to characterize rock cores? T 2 is considered the industry standard metric, the reasons for this being (i) that T 2 represents a more local measurement of the pore-size distribution because it is averaged over a shorter timescale measurement than T 1 and (ii) that only T 2 (and, to a lesser degree, D) measurements are carried out with MR well logging tools. The longer timescale of the T 1 measurement means that the acquired data will be averaged over more than a single pore. So, given the limitation that T 2 varies significantly with field strength, why would we wish to make measurements in higher field magnets? There are two reasons: (i) improved SNR at higher field strength means faster data acquisitions, important for studying time-variant systems, and (ii) higher field strengths offer the opportunity for spectral resolution, useful for simultaneous quantitative identification of oil and water fractions. Having answered these questions, we can now pose the following question: is it possible to attain meaningful T 2 distributions at high fields on relatively MR 'unfriendly' materials such as sandstones in which the pores might be lined with iron-bearing clays? We note that a prerequisite of any such T 2 measurement is that the derived relaxation time distribution must be consistent with the poresize distribution of that material; a method of achieving this was described in section 2.1. Here, we highlight another element of that work which employs Bentheimer as a typical sandstone for the purposes of demonstrating the complexities of acquiring accurate T 2 distributions in highfield spectrometers (Mitchell et al 2010a), where the magnetic susceptibility contrast between the rock and brine is χ = 2 × 10 −6 (in SI units). In this study, T 1 and T 2 data were acquired over a range of field strengths from B 0 = 0.042 to 9.4 T. At low magnetic fields (B 0 < 1 T; ST regime), the maximum effective internal gradient strength scales approximately as B 3/2 0 , although at high magnetic fields (B 0 > 1 T; the LOC regime), this trend is limited by the minimum pore size (Mitchell et al 2010b). Therefore, the maximum effective gradients in the LOC regime may not be as large as predicted by theory based on the ST regime (see section 2.1). Bentheimer has a monomodal pore-size distribution (determined by mercury intrusion porosimetry) and this was reflected in T 1 distributions acquired across the range of field strengths considered. The log-mean T 1 increased with increasing field strength as predicted theoretically. However, it was found that the equivalent (effective) T 2,eff distributions became more complicated (i.e. multimodal) as the field strength increased above approximately B 0 = 1 T. This multimodality was accompanied by a reduction in log-mean T 2,eff , contrary to theoretical predictions for T 2 . This variation in T 2,eff was observed dramatically in 2D T 2,eff -T 2,eff exchange measurements (see figure 15). As the magnetic field strength increased, more onand off-diagonal peaks were observed in the exchange measurement. The off-diagonal peaks were attributed to molecular diffusion of water between regions of differing internal gradient strengths within individual pores. The average exchange rate corresponded to d ∼ 60 µm: the mean pore size based on mercury intrusion and electron microscopy measurements. It was therefore proposed that at high-magnetic field strengths and/or in the presence of strong internal gradients, T 2 -T 2 exchange measurements provide an alternative measure of mean pore size.
The use of MR propagator measurements to characterize flow in rock cores is well established. However, to make these measurements quantitative as opposed to qualitative remains challenging. Scheven et al (2005) measured pre-asymptotic Stokes flow through (i) a bead pack with a bead diameter d = 100 ± 20 µm, (ii) a Bentheimer sandstone and (iii) a Portland carbonate rock core, for Peclet numbers in the range Pe = 20-80 and Reynolds numbers of Re < 0.1. Probability distributions of molecular displacements P(ζ ) were analyzed by calculation of the mean displacement ζ , the rms width σ and adjusted skewness γ using the self-consistent cumulant analysis designed to minimize systematic errors, as described in section 2.2. Figure 16 shows a selection of the propagators obtained, plotted in the nondimensionalized form P(ζ ) for ζ > 0, where ζ 0 is the expected mean displacement along the direction of superficial flow. Propagators for five observation times in the range = 106-2000 ms are shown. It is seen clearly that the shape of the propagators for flow within the bead pack evolves by mechanical mixing from an initially skewed pre-asymptotic shape (for small mean displacements) to one resembling a Gaussian centred about the mean displacement (for large displacements). The decreasing width of the distribution as mean displacement increases is consistent with mechanical mixing and it is also clear that there is no evidence of stagnant fluid within the pore structure, as expected for this type of regular packing. The  (Scheven et al (2005)). sandstone and carbonate show very different behaviour. The distributions remain skewed at all values of but do show evidence of mechanical dispersion. Both exhibit pronounced peaks around zero displacement: clear evidence of stagnant zones. The flow field in the carbonate is much more heterogeneous than that of the sandstone. The stagnation peaks are larger and faster flows exist within the carbonate. Despite the use of propagators for 'fingerprinting' the transport characteristics of a particular fluid within the rock, as noted in Scheven et al (2005), these displacement distributions are distorted from their true form. Particular contributions to the distortion arise from both slow and fast moving spins (see section 2.2). Mitchell et al (2008b) have considered analysis of flow propagators obtained for Bentheimer sandstone and Portland carbonate using both the analytic correction of Scheven et al (2005) and their own empirical approach using a combination of ζ -T 2 and T 1 -T 2 correlations. The displacement l corresponding to equivalence in the first and second moments (and hence providing an indication of the dispersivity) was determined. In the uncorrected data, l = 434 and 169 µm for the carbonate and sandstone, respectively. Following the analytic correction of Scheven et al (known to overestimate the fraction of stagnant spins), the displacements increased to l = 1054 and 209 µm. However, the empirical relaxation correction provided displacements of l * = 758 and 199 µm for the carbonate and sandstone, respectively; so l < l * < l as expected. From this study, we see that if accurate displacement distances are to be determined, the empirical relaxation correction should be applied. However, these experiments are time consuming and so, failing this, an improved displacement distance can be obtained from an average of the uncorrected and analytically corrected cases.
Bearing in mind the techniques we can now employ to improve our measurements of T 2 , T 1 -T 2 , T 2 -T 2 and propagator displacement distributions, it is useful to remember that the addition of spectral resolution into T 1 -T 2 correlations has also been demonstrated recently, which offers opportunities for estimating the surface wettability in rock cores (Chandrasekera et al 2008, Mitchell et al 2009b. Figure 17 shows the results of a T 2 -T 1 -δ correlation for a 50 : 50 mix of water and dodecane in two carbonates-one that is water wetting and a second that has been treated to become partially oil wetting. The water-wet case is shown in figures 17(a) and (b); the chemically resolved water signal is shown in (a) and is seen to be shifted to a higher T 1 /T 2 ratio as expected due to its surface interaction. The oil signal (b) remains on the T 1 = T 2 line (consistent with bulk dodecane), evidence for the oil not interacting with the pore surface. In the partially oil-wet case (figures 17(c) and (d)), both the water and oil show evidence of interaction with the surface.

Processing-structure-function relationships: cement and gypsum.
MR is a particularly powerful technique in studying aspects of manufacturing processes, the influence of different process parameters on the microstructure of the resultant materials and the effect of that microstructure on the in-use properties of the material. We have taken recent studies of cement and gypsum to illustrate emerging applications of MR in these areas.
2D T 1 -T 2 relaxation correlations have been used to study the evolution of hydrating white (architectural) cement pastes (McDonald et al 2005). Two types of pore structure form in cement during hydration: small calcium silicate hydrate (CSH) gel pores (d = 1-30 nm) evolve as a consequence of the crystallization process, while large capillary pores (d ∼ 100 nm-10 µm) form due to chemical shrinkage. The T 1 -T 2 correlation plots showed a series of distinct relaxation time peaks all exhibiting the ratio T 1 /T 2 = 4 at a 1 H Larmor frequency of f 0 = 20 MHz. These peaks were attributed to two sizes of gel pores (T 1 = 200 µs and 2 ms) and capillary pores (T 1 > 10 ms). An additional off-diagonal peak, attributed to exchange between the two sizes of gel pore, was observed with T 1 = 10 ms and T 2 = 50 µs. This exchange was probed further with T 2 -T 2 correlations from which a diffusive exchange rate of k = 5 ms −1 was obtained for water moving between the gel pores (Monteilhet et al 2006). T 1 -T 2 correlations of hydrating cements have also been obtained using a unilateral magnet (McDonald et al 2007a) for both white cement and Portland (grey) cement that are used widely in the construction industry (McDonald et al 2007b).
Similar approaches have been taken to study gypsum plaster pastes. Here, we seek to understand the relationship between the manufacturing process and the physical properties (i.e. strength and water resistance) of the resulting plasterboard. A particular area of interest has been monitoring the kinetics of hydration and the evolution of pore structure during this process; this is more challenging than in cements because the timescales of hydration in plaster paste are significantly shorter (in the order of minutes compared to days). T 2 relaxation time measurements have been used to monitor the evolving pore structure in β plaster (Jaffel et al 2006a). Jaffel and co-workers observed two distinct water populations, and the variations in the relative water fraction of these populations was monitored as a function of hydration time for a range of water-to-plaster ratios (w/p). The pore structures determined from the MR data were in agreement with SEM images of the set products, and were found to be related to the macroscopic mechanical properties as determined by ultrasound velocity measurements (Jaffel et al 2006b). Elsewhere, it has been shown that the hydration and deterioration of wet plaster can be followed using T 2 relaxation measurements on a low-field, portable MR device (Boguszynska et al 2005). Further work has compared the hydration kinetics of α and β plaster pastes using 1D NMR spin density profiles, and the evolution of the pore size has been studied by T 2 relaxometry (Song et al 2009). In this work, the relaxometry data show clear evidence of the formation of large pores within the β plaster, which are not present in the α plaster-a consequence of different treatments used during the production of the dry powder. The theme of studying manufacturing-structure-property relationships is being explored further with aerated gypsum products. Aerated gypsum is produced by incorporating pre-generated foam directly into the wet gypsum slurry, which leads to the formation of air voids within the plasterboard structure. In addition to characterizing the pore structure using T 2 relaxometry, water imbibition into the set plaster has been studied. This revealed the presence of two water penetrating fronts moving through the material. This wetting behaviour can, of course, be related back to the microstructure characterization. A detailed consideration of the determination of sorptivity in these materials using MR, and comparison with sorptivity data obtained using conventional gravimetric techniques, has been presented by Song et al (2011). MR is shown to provide far greater insight into the details of the microporous structure within the gypsum. In particular, the presence of two wetting fronts due to the existence of a bimodal pore-size distribution is evident in the MRI data but could be overlooked easily in the gravimetric measurements. The pore-size distribution and connectivity of pores of different size were shown to have a significant impact on water transport within the plasterboard, and therefore to influence the in situ performance of the finished products.
Finally, rapid 2D T 1 -T 2 correlations (see section 2.1) have been employed in a novel study of additives included in gypsum plaster pastes to modify (accelerate or retard) the hydration processes. These additives are presumed to function by modifying the interaction of water with the evolving gypsum crystal structures and hence influence the resultant pore network. The T 1 /T 2 ratio was therefore utilized to provide an indication of both the number and strength of interactions between water molecules and solid gypsum particles (Song et al 2010). Relative to the T 1 /T 2 ratio of unmodified gypsum paste, ground gypsum mineral (GGM) increased the T 1 /T 2 ratio, whilst citric acid decreased the T 1 /T 2 ratio. This was interpreted as an increase in the water-gypsum interaction in the presence of GGM, and a decrease in the presence of citric acid, consistent with the current understanding of the behaviour of these additives as an accelerator and retardant, respectively.

Catalysts.
Most MR imaging, diffusion and surface interaction studies have focused on supported metal catalysts or the support materials themselves (e.g. porous silica, alumina or titania) as opposed to microporous crystalline materials, such as zeolites and other molecular sieve catalysts. In the context of this review, we note that PFG NMR has been used to study diffusion in zeolites (Datema et al 1991, Feldhoff et al 2009, although care must be taken if quantitative data are to be acquired (Kärger and Vasenkov 2005). 2 H MR has also been used to characterize intra-cage and inter-cage transport of aromatics in Y-zeolites (Sousa Gonçalves et al 1995).
The majority of early MRI studies specific to catalysis addressed the heterogeneity in structure and transport within oxide catalyst pellets. In-plane spatial resolution achieved in these studies was of the order of 30-50 µm for an image slice thickness of 0.3 mm and the pellets themselves were of typical dimension 1-5 mm. Spatially resolved maps of spin density, T 1 and D of water within a porous catalyst support pellet (Hollewand and Gladden 1993, 1995a, 1995b were reported with in-plane spatial resolution of 45 µm × 45 µm and an image slice thickness of 0.3 mm. Estimates of overall pellet voidage obtained from spin density maps agreed to within 5% of that obtained by gravimetric analysis. Even though the coarse resolution T 1 maps provide only an estimate of mean pore size within a given image pixel, these data give us insights additional to those from bulk pore-size distributions as obtained, for example, by mercury porosimetry or nitrogen adsorption. In particular, by using MR, it became possible to probe the spatial heterogeneity in porosity within a catalyst pellet that had been introduced during the manufacturing process. There have been a number of studies exploiting this type of MR measurement. For example, water spin density imaging studies have been used to explore the 3D structure of activated alumina spheres-the spheres were observed to exhibit a uniform ball structure, comprising spherically layered sub-structures and voids (Timonen et al 1995). Simple spin-density imaging has also been used to characterize the tortuosity of catalyst pellets manufactured by different processes . This was achieved by impregnating the catalyst initially with deuterated water (invisible to a 1 H MR 31 experiment) and then immersing it in normal protonated water. The diffusive exchange of 1 H 2 O with 2 H 2 O within the pellet was then followed by 1 H imaging. The effective water diffusivity and hence catalyst tortuosity were obtained by fitting the time-resolved 1 H 2 O concentration profiles to a standard diffusion equation. Measurements of this type are straightforward to perform and give immediate insight into transport anisotropies within the catalyst resulting from the manufacturing process. Later work demonstrated that it is possible to characterize the mesoscale porosity of a given batch of catalyst pellets by a fractal dimension . As mentioned earlier, numerical modelling strategies employing simulation lattices constructed to have the same fractal dimension as those obtained from the MR images have shown some success in predicting mass transfer characteristics in porous catalyst supports (Rigby and Gladden 1998). This toolkit of MR measurements, namely 1 H spin density imaging with the ability to exploit relaxation contrast, has also been used to study aspects of catalyst preparation and drying behaviour of catalysts with reference to catalyst manufacture , Espinosa-Alonso et al 2009.
MR has also been applied to pellet-scale imaging of the effects of reaction occurring within catalyst pellets. Koptyug et al (2002) have reported 1 H images of liquid distribution during α-methylstyrene evaporation accompanied by its vapour-phase hydrogenation within a cylindrical Pt/γ -Al 2 O 3 catalyst pellet (diameter and height 4.7 mm). Imaging and relaxometry methods have also been used to probe how hydrocarbon (coke) residues are deposited within the catalyst, thereby modifying its surface chemistry and pore structure and therefore influencing transport and reaction processes. A number of different approaches have been reported to monitor the spatial distribution of coke and the qualitative evolution of the carbon-to-hydrogen ratio characterizing that coke within porous alumina pellets and zeolites (Cheah et al 1994, Bonardet et al 1999, Bär et al 2002, as well as the effect of coke deposition on the self diffusion of hydrocarbons within industrial catalysts (Wood and Gladden 2003).
A more recent extension of the ability to acquire spatially resolved diffusion maps (Hollewand and Gladden 1993) is shown in figure 18. Here, a chemically selective diffusion mapping sequence was used to identify simultaneously the mobility of 2-butanone and 2-butanol in a single spherical Ru/SiO 2 catalyst pellet (Lim et al 2011). The apparent diffusivity of each liquid inside the pellet was reduced by a factor of 2 approximately, compared to the bulk liquid diffusivity. The motivation for this work is to develop the measurement capability to probe gradients in concentration and molecular mobility within individual catalyst pellets so that we can characterize and understand the origin of mass transfer resistance in catalyst pellets during reaction.
Finally, two examples are now given of how MR techniques are being used to probe molecular dynamics on the internal surface of the catalyst. First, T 1 -T 2 relaxometry techniques-described above in reference to surface wetting in rock core analysis-are now being applied to catalytic materials to study the behaviour of liquid phase reactants and products within the pore space. As demonstrated earlier, the relative increase in the average T 1 /T 2 ratio when a molecular species adsorbs on a surface gives an indication of the strength of that surface interaction. In catalysis, we are interested in the same wetting phenomenon: in a multicomponent liquid system, which species dominates the catalyst surface? This is of great importance in understanding and optimizing catalysed reactions: do the reactants adsorb and react or do they get blocked from accessing the surface by product or solvent species? Weber et al (2009)  hydrogenation of 2-butanone (namely 2-butanone, 2-propanol and water) with Ru/SiO 2 and Pd/Al 2 O 3 . Secondly, PFG NMR has been used to study diffusion inside porous catalysts. Weber et al (2010) used an APGSTE sequence with the gradient amplitude varying between g = 0 and 12 T m −1 and observation times of = 10,75, 150 and 250 ms, to probe diffusion in liquidsaturated catalyst pellets. For the first time, evidence was seen for the existence of two diffusion regimes, attributed to bulk liquid diffusion in the catalyst pores and restricted diffusion of molecules adsorbed on the pore surfaces. By fitting a model for two-site exchange to these data the values obtained for bulk pore and pore-surface diffusion for 1-octene within a 1 wt.% Pd/Al 2 O 3 catalyst were found to be 1.3 × 10 −9 and 1.7 × 10 −11 m 2 s −1 , respectively, with the estimated mean residence time of a molecule on the pore surface being τ s = 150 ms. Data of this nature were also recorded for pure 2-propanol and pure water within porous alumina pellets, and the results were compared with molecular dynamics simulations (Youngs et al 2009). Trends in the MR and molecular dynamics results were consistent and a direct comparison suggested that the MR measurement of surface diffusion probes two molecular layers adjacent to the catalyst surface. This result is consistent with the chemical engineers' 'rule-of-thumb' regarding the extent of the surface-influenced layer, and this comparison of MR experiment with molecular dynamics simulation provided a useful confirmation of this assumption.

Macroporous systems
3.2.1. Fixed beds. The MR techniques used to investigate macroporous systems focus on flow imaging, propagator measurements and, increasingly, chemical mapping (i.e. spatially resolved chemical composition). Many of the systems of interest are also characterized by porosity over a hierarchy of length scales with a microporous network interconnected with mesoporous and/or macroporous networks. Hence, we are not interested solely in the flow through the macropores but are also concerned with how this flow communicates with networks of smaller pores. A more detailed review of MR applications in catalysis and reaction engineering is found elsewhere . This brief overview of recent developments will address fixed-beds of packing elements of direct relevance to heterogeneous catalytic processes. However, the approaches taken have been readily spun out in many other applications, including bioseparations (Holland et al 2004), biofilms (Graf von der Schulenburg et al 2007a, 2008, Pintelon et al 2010 and dissolution testing of controlled release pharmaceutical systems (Shiko et al 2010).
High-resolution MRI studies of single-phase fluid flow in packed beds were first reported in the mid-1990s. These early studies were performed on beds with column-to-particle diameter ratios in the range 10-20, typical of narrow fixed-bed reactors and used non-porous packings, e.g. solid glass spheres (Sederman et al 1997(Sederman et al , 1998. From spatially resolved measurements of flow velocity in the x-, yand z-directions, it is possible to calculate maps of shear stress within the fluid (Sederman and Gladden 2001a). These maps show precisely where high liquid shear stress exists within the bed; such regions will be associated with particle erosion but are not likely to be regions in which significant fines deposition will occur. These regions are also associated with high fluid velocities and inertial effects that increasingly influence the flow profile (Johns et al 2000). In contrast, in regions of low shear stress, little particle erosion will occur but we will expect these regions to be particularly susceptible to fines deposition. The type of information obtained from these images can be used both at an academic level and by catalyst manufacturers interested in investigating how subtle changes in pellet size and shape can influence the spatial distribution of porosity and hence hydrodynamics within the reactor (e.g. Götz et al 2002). Suekane et al (2003) reported MR velocity imaging studies of flow up to Re = 204.74 in a narrow packed bed and investigated the relative importance of inertial to viscous forces as a function of Re using the concept of a local Reynolds number. MR radial velocity profiles of single-phase flow in narrow fixed beds with column-to-particle ratio <10 have also been combined with traditional reaction engineering characterization of the same beds, namely residence time distribution analysis (Tang et al 2004). In particular, the influence of column-to-particle diameter and the Reynolds number (based on particle diameter) on dispersion and the degree of non-uniformity of the velocity profile were studied. Unsteadystate flows in fixed beds have also been imaged; the GERVAIS pulse sequence has been applied to image unsteady flow in a fixed bed with a column-to-particle ratio of 2 (Sains et al 2005).
Two-phase gas-liquid flow has also been studied with particular attention paid to using MRI for gaining insight into the hydrodynamic transition from a steady-state gas-liquid distribution at low gas and liquid velocities (trickle flow) to the unsteady state gas-liquid distribution referred to as pulsing flow. However, let us consider first the trickle flow regime. In general, when studying trickle flow, high spatial resolution images (i.e. of the order of 100-200 µm in-plane resolution) are used, typically, to obtain (i) accurate measurements of liquid and gas holdup within the bed where the liquid holdup is in quantitative agreement with that measured gravimetrically and (ii) direct measurements of liquid-solid contact (i.e. catalyst wetting) Gladden 2001b, Gladden et al 2003a). More recently, we have extended our imaging capability to obtain both liquid and gas velocity maps in fixed beds. Sankey et al (2009) reported 3D velocity images of liquid flow in the superficial flow direction in a 27 mm diameter bed packed with 5 mm glass spheres. Flow images were acquired at an isotropic spatial resolution of 266 µm. This work also reported the first images of both Figure 19. Gas and liquid velocity maps of SF 6 (red/yellow) and water (blue/green) during trickle flow. The gas and liquid superficial velocities were 8.7 and 2.3 mm s −1 , respectively. The bed was operated at a pressure of 4.7 bara. Reprinted with permission from Elsevier (Sankey et al (2009)).
gas and liquid velocity within a packing of spheres; the liquid and gas were water and SF 6 , respectively. The gas velocity map was acquired using 19 F observation. An example of the data obtained is given in figure 19. The motivation for acquiring gas and liquid velocity maps within the same bed is to gain insight into the closure relationships needed in the implementation of computational fluid dynamics codes. For example, we could ask: what is the relationship between the gas and liquid across a gas-liquid interface? At the moment, such relationships are often derived empirically. Despite the success of acquiring the data shown in figure 19, there remains a problem with regard to extracting the information required to answer this question. This is because the spatial resolution achievable in the gas velocity map is far lower than that of the liquid velocity map as a result of the much lower gas density. In the data shown in figure 19, the liquid map is taken for a 1 mm slice thickness with in-plane spatial resolution of 177 µm × 177 µm. In contrast, the gas image is acquired from a 2 mm slice thickness at an in-plane spatial resolution of 708 µm × 708 µm. Therefore, we cannot avoid mis-registration of the gas-liquid interface that is exactly the region we want to study. Overcoming this problem was a major motivation for implementing CS methodologies into our data acquisitions. Holland et al (2010a) have demonstrated recently the application of CS to single-phase gas and liquid flows. Using sampling strategies in which ∼30% of the k-space is acquired has enabled both gas and liquid flow images in the system described above to be acquired at the same slice thickness and spatial resolution: namely a slice thickness of 1.5 mm and an in-plane spatial resolution of 230 µm × 230 µm. Work is ongoing to analyse these data.
Another area where MR has been shown to be particularly useful in reaction engineering is in identifying the mechanism of the hydrodynamic transition from trickle to pulsing flow. 2D and 3D fast low-angle shot (FLASH or SNAPSHOT) pulse sequences (Haase et al 1986) were used to image the temporal stability of liquid distribution within a fixed bed under conditions of gas-liquid two-phase flow (Lim et al 2004, Anadon et al 2006. The transition was seen to develop through the formation of increasingly large numbers of local pulsatile events, which eventually formed a spanning cluster within the bed. At this point, the bed is described as 'pulsing'. Much of the work focused on a fixed bed of length 700 mm and Time-dependent behaviour of local pulsing events. 3D standard deviation maps combined with a RARE image of the bed structure are shown for a bed height of 34 mm. All maps were calculated from 3D MRI data sets acquired with the bed operating under the same conditions with gas and liquid velocities being 300 and 6.1 mm s −1 , respectively. The data for each map were acquired in less than 2.24 s and the time interval between the data acquisition for each of the standard deviation maps was 2 s. Regions A and B identify volumes within the bed at which local pulses are observed; these local pulses 'switch on' and 'switch off' as a function of time. There also appear to be regions in the bed, such as that identified as C, in which local pulses are not observed, suggesting that the local structure of the bed plays a role in determining the spatial location of the local pulsing events. High values of standard deviation (indicating unstable liquid content) are identified as red, and low values (indicating stable liquid content) are identified as blue/green. Reprinted with permission from John Wiley and Sons (Anadon et al (2006)).
of diameter 43 mm. The need to identify time-varying liquid content within the bed meant that temporal resolution was far more important than spatial resolution in these experiments. To study instabilities in gas-liquid distribution in 2D, 540 images were acquired. Each image had a data acquisition time of 20 ms, and images were acquired in immediate succession. For the case of the 3D data, six series of eight consecutive 3D images were acquired for any given combination of gas and liquid flow rates. Superficial gas and liquid velocities in the range 25-300 and 0.9-13.0 mm s −1 , respectively, were used. Each 3D image was acquired in 280 ms (Anadon et al 2006). The variation in liquid content within each pixel (voxel) over the complete time series of images acquired was achieved by calculating the standard deviation in pixel (voxel) intensity on a pixel-by-pixel basis and then displaying the spatially resolved map of standard deviations in 2D or 3D. Regions characterized by low values of standard deviation identify stable gas-liquid distributions, while values of high standard deviation identify spatial locations in which the liquid content (i.e. gas-liquid distribution) is fluctuating with time. Figure 20 shows standard deviation maps calculated from the 3D images of time-varying liquid distribution. A cross-section through the bed reveals the evolution of the local pulses in a direction parallel to the direction of superficial flow. Further MR experiments suggest that fluctuations in the thickness of liquid films on the surface of the packing elements increased with increasing liquid velocity prior to the formation of local pulsing events (Lim et al 2004. The opportunity to translate the understanding of this hydrodynamic transition gained by MR such that cheap, robust sensors can be used on industrial-scale chemical reactors has also been explored (Anadon et al 2008).
Chemical mapping in a fixed-bed reactor was first demonstrated by Yuen et al (2002) using the liquid phase esterification of methanol and acetic acid catalysed within a fixed bed of H + -ion exchange resin (particle size 600-850 µm) using volume selective spectroscopy (Kimmich and Hoepfel 1987). The ability to estimate mass transfer coefficients in situ and to explore correlations between local conversion and hydrodynamics has also been demonstrated in application to this esterification reaction (Gladden et al 2003b, Yuen et al 2003. However, to study the hydrocarbon conversions more typical of the chemical industry, 13 C observation is needed so as to benefit from the increased chemical shift range such that reactants, intermediates and products can be identified. Thus, 13 C DEPT MRI has been employed and was demonstrated initially in the study of a competing etherification and hydration reaction of 2-methyl-2-butene (Akpa et al 2005). Following this work, 13 C DEPT MRI was applied to map isomerization and full conversion of 1-octene during hydrogenation over a Pd/Al 2 O 3 catalyst . Quantitative interpretation of these data is challenging because of the influences on the spectral line shape from (i) species existing in both the intra-and inter-particle space as well as (ii) individual peaks within the spectrum of a given molecular species being attenuated to differing extents due to their different interactions with the catalyst surface. This has led to considerable effort to implement a partial least squares (PLS) chemometrics approach to analyse objectively the spectra we obtain. Data can be acquired sufficiently fast such that the spatially resolved composition in the bed can be followed in unsteady-state systems. In figure 21, the mole fractions of 1-octene, 2-octene and n-octane along the length of the reactor are shown at two times: 22.5 and 180 min after the introduction of hydrogen. The time associated with a given data set is the time at the half-way point through the total data acquisition time. From these data, the conversion and selectivity to 2-octene and n-octane can be obtained as a function of both time and axial position along the bed. Measurements of the individual components of granular temperature for nudicaule seeds: (a) T x (b) T y (c) T z . Note that the scale has been increased for the image in (c). The superficial gas velocity was 0.26 m s −1 . The spatial resolution of the images was 2.50 mm (z) × 1.45 mm (x) and the slice was located in the centre of the bed with a thickness (y-direction) of 5 mm.

Fluidized beds.
Imaging the motion and distribution of solids directly in a gas-solid fluidized bed is another area of increasing activity (Savelsberg et al 2002, Fennell et al 2005, Harms et al 2006. Different approaches may be taken to acquire signal from the solid phase. The most common approach is to use edible solid particles, such as poppy seeds, in which the natural oil content provides a strong 1 H MR signal. Alternatively, a porous solid can be impregnated with a liquid (containing 1 H species typically) to provide the MR signal. A number of studies have now been published that explore the use of MR in this field. For example, 1D FLASH imaging has been used to measure simultaneously both the rise velocities and coalescence of bubbles, and the dynamics of the solid phase, in gas-solid twophase flow. Experiments were performed with a bed of particles of diameter 0.5 mm contained within a column of internal diameter 50 mm. Gas velocities in the range 0.18-0.54 m s −1 were studied. Data were of sufficient temporal and spatial resolution that bubble size, and the evolution of bubble size and velocity following coalescence, were determined. Spatially resolved measurements of granular temperature have also been reported ; typical results are shown in figure 22. Granular temperature arises from an analogy between the motion of individual particles in granular systems and the motion of molecules in kinetic theory. Collisions between individual particles will cause random fluctuations in the particle velocity, analogous to the thermal fluctuations in molecules; the magnitude of these velocity fluctuations is termed 'granular temperature'. Spatially resolved velocity maps of solids motion allows direct calculation of local granular temperature within the fluidized bed.
Continuing the philosophy underpinning the fixed-bed work, it has also been possible to explore how MR measurements relate to and complement the results from other sensor and tomographic methods. This is of particular importance for work on gas-fluidized beds because these process units will operate at much larger size scales than it is possible to study within the confines of a superconducting magnet. We note that pilot units for fixed-bed processes are at the size scale studied using commercial MRI hardware. Examples of comparisons between MR and other methods include the study of oscillations in gas-fluidized beds, in which ultrafast MRI data have been compared with pressure sensor measurements and used to assess the predictive power of existing literature correlations for pressure fluctuations within the bed (Müller et al 2007). In a more recent example, MR data have been compared with electrical capacitance volume tomography (ECVT) data . Snapshot MR techniques were used to acquire 2D data at a temporal resolution of 26 ms with in-plane resolution of 1.9 mm × 1.9 mm in the transverse plane and 1.9 m × 3.8 mm in the axial plane. The ECVT data were acquired at a temporal resolution of 12.5 ms with a nominal spatial resolution of 2.5 mm × 2.5 mm × 4.5 mm. The two techniques produced quantitatively comparable bubble frequencies and time-averaged measurements of voidage. The definition of the gas bubble and the wake structure were resolved more easily in the MR data.
The wealth of data that can now be obtained on fluidized beds-i.e. solids distribution, bubble velocity, coalescence and formation of larger bubbles and slugs, velocity maps and granular temperature-provide a much needed resource for the development and validation of numerical codes. In the case of granular flows, discrete element modelling is the approach of choice. Recent work  has demonstrated good agreement between the predictions of a discrete element model simulation with a map of granular temperature obtained using MRI. This work is a good example of how MRI data can be used in combination with code development to help to identify the correct implementation of the code to give agreement with the image data.

Conclusions
Over the past 20 years, the application of MR techniques to the study of chemical engineering processes and process-structure-function properties of materials has expanded significantly. One of the challenges of applying MR to porous media and related systems encountered in chemical engineering is the existence of rapid T 2 and T * 2 relaxation times and the existence of internal field gradients caused by variations in magnetic susceptibility within the system of interest. In applying MR to porous systems in chemical engineering, the acquisition of quantitative data that can be used for comparison with numerical or theoretical predictions or as input to such simulations and theories is of utmost importance. Strategies for improving the quantitative nature of measurements of MR propagators is a good example of how measurement protocols to remove relaxation weighting in complex porous structures has been advanced in recent years. In conclusion to this review, we ask: what do we need, and therefore expect, from future developments? Our focus is on two particular areas: (i) spatially resolved chemical resolution at faster data acquisition rates and how we can exploit polarization enhancement techniques without losing the quantitative nature of the MR signal, and (ii) faster data acquisitions through sparse k-space sampling that will not only enable us to probe systems fluctuating on shorter timescales but also allow us to perform greater data averaging, thereby increasing the SNR available in a given data acquisition time. In this regard, the areas that we are interested in developing are CS and Bayesian MR. We also take this opportunity to highlight the importance of designing experiments which contribute relevant information to real chemical engineering problems. For example, experience shows that it is important to understand transport in the pores of rock cores when optimizing oil recovery processes, even though in the field, macroscale heterogeneity of rock formations will clearly influence the process performance. Considering the macroporous systems reviewed here, we note that fixed beds have been a major topic of our studies since industrial scale-up is often performed based on beds of the size we can accommodate readily within our magnets. In contrast, we expect the information that we can access on centimetre-scale gas-fluidized beds to give useful information about local bubble dynamics and to provide validation of numerical codes, but we do not expect to obtain information relevant to the macroscale hydrodynamics typical of commercial beds constructed at metre-scale dimensions. In such applications, we will rely on the numerical codes-validated by MR-to be used in scale up studies. Overall, this review has highlighted the significant steps taken so far in MR applications to chemical engineering. We have presented the foundations that have been laid for the next generation of MR pulse sequence development and given some insight into the future directions of industrial applications of MR. Despite all of the successes, there are still many more exciting applications to be explored by the ever-expanding toolbox of MR measurements, and increasing combinations of MR techniques with numerical simulations and other tomographic methods. The versatility of MR to probe transport and interactions across the hierarchy of length scales relevant to industrial process applications will ensure the continued importance and progress of these techniques as the chemical engineering solution of choice for studying many porous systems.