Real-time probing of granular dynamics with magnetic resonance

See allHide authors and affiliations

Science Advances  15 Sep 2017:
Vol. 3, no. 9, e1701879
DOI: 10.1126/sciadv.1701879


Granular dynamics govern earthquakes, avalanches, and landslides and are of fundamental importance in a variety of industries ranging from energy to pharmaceuticals to agriculture. Nonetheless, our understanding of the underlying physics is poor because we lack spatially and temporally resolved experimental measurements of internal grain motion. We introduce a magnetic resonance imaging methodology that provides internal granular velocity measurements that are four orders of magnitude faster compared to previous work. The technique is based on a concerted interplay of scan acceleration and materials engineering. Real-time probing of granular dynamics is explored in single- and two-phase systems, providing fresh insight into bubble dynamics and the propagation of shock waves upon impact of an intruder. We anticipate that the methodology outlined here will enable advances in understanding the propagation of seismic activity, the jamming transition, or the rheology and dynamics of dense suspensions.


Traditionally, the rich and industrially relevant physics of dynamic granular systems have been investigated by inserting probes into three-dimensional (3D) flow systems. Because of the nature of probes to interfere with the flow dynamics, noninvasive techniques for measuring granular motion have been developed. Optical techniques can image particle distribution and dynamics but are limited to measurements in pseudo 2D systems and close to the walls of 3D systems. More recently, advanced tomographic techniques, such as x-ray imaging (1) and tomography (2), positron emission particle tracking (3), electrical capacitance tomography (4), and magnetic resonance imaging (MRI) (516), have been used to characterize granular dynamics throughout the volume of full 3D systems (17). Most techniques have been limited to measurements of the distribution of particles (solids fraction), but MRI is far more versatile. MRI has allowed for quantitative and spatially resolved measurements of particle distribution and velocity (58), acceleration (10) and diffusion (11), force chain networks (12), temperature (13), the progress of chemical reactions (14, 18, 19), and the dynamics of interstitial gas (15), thus making it a promising technique to study the complex dynamics of single- and multiphase granular systems (512, 15, 16, 20, 21).

Two key disadvantages of MRI have been small system sizes and low spatiotemporal resolution. Granular systems studied using MRI have been limited in size such that flow features are strongly affected by walls and dynamics are simplified as compared to industrial-scale systems (5, 6, 8). Granular systems involve complex 3D dynamics, with flow features changing significantly over time scales down to milliseconds and length scales comparable to the particle diameter (22). Limits in spatiotemporal resolution have only allowed for time-averaged (>102 s) 2D images of particle velocity (58) or rapid 1D measurements (23). Resolution limitations largely stem from the fact that to obtain images with sufficient signal-to-noise ratio (SNR), several excitations of nuclei per image were required, with significant delays between excitations.

Here, we present an ultrafast single-shot MRI methodology, which produces images of particle distribution every 7 ms in systems large enough to allow for complex, disordered flow using only one excitation pulse (single-shot MRI). We extend this technique further to enable the imaging of particle velocity with a temporal resolution of 22.5 ms, effectively advancing MRI velocimetry of granular dynamics from temporally averaged images to instantaneous snapshots. Furthermore, we present a new contrast mechanism to detect propagation of waves of coherent particle motion in granular media.

This single-shot methodology was enabled by reducing the image acquisition time while simultaneously increasing the signal from the particles during acquisition. First, array detection of the MR was applied by placing 16 radio frequency (RF) receiver coils around the region of interest. The spatially varying sensitivities of the receiver coils allow for a reduction in the sampling density of k-space below the Nyquist limit (24), effectively decreasing the acquisition time (Tacq) compared to single-channel acquisition (25). Second, a partial Fourier sampling scheme combined with a homodyne phase–constrained image reconstruction (26) was used to further reduce Tacq by acquiring only slightly more than half of k-space (lower green part of Fig. 1D). Third, inspired by previous granular MRI studies using pharmaceutical pills filled with vitamin E (27, 28), spherical particles composed of an agar shell filled with oil were engineered (particle diameter dp = 1.02 ± 0.12 mm). The advantageous MRI characteristics of these particles increased both the initial signal arising from the particles and the signal lifetime (T2*) as compared to conventionally used agricultural seeds. The interplay between these three advances (TacqT2*) enables the time-efficient acquisition of a full image from a single excitation with sufficient SNR in granular systems.

Fig. 1 Enabling single-shot MRI of granular dynamics.

(A) Concurrent detection of MR using an array of 16 RF coils placed closely around the region of interest. The varying spatial sensitivities of the coil allow for a reduced k-space sampling density below the Nyquist limit (see Materials and Methods). (B) Single-shot echo-planar imaging (EPI) readouts are applied to sample k-space time-efficiently within one spin coherence. Asymmetric partial Fourier sampling reduces the readout duration further. The green shaded area of k-space is acquired, whereas the gray shaded area is estimated in the algebraic image reconstruction (see Materials and Methods). (C) Photograph of tailored oil-filled agar particles (particle diameter dp = 1.02 ± 0.12 mm) with enhanced MR signal intensity and lifetime when compared to agricultural seeds provides sufficient signal throughout the single-shot readout duration of 4.6 ms. (D) k-space representation of the single-shot EPI readout trajectory (bold black dots). The light dashed line represents a conventional EPI sequence sampled according to the Nyquist limit. The zoomed area shows the applied reduction of sampling density. The pulse sequence properties for this readout trajectory are detailed in fig. S1A. (E) Algebraic image reconstruction exploiting varying spatial sensitivities of the 16 receive channels to correct for undersampled, asymmetric k-space information. (F) Reconstructed image obtained by single-shot MRI. The temporal and spatial resolution (Δx × Δy × Δz) of the acquisition was 7 ms and 3 mm × 5 mm × 10 mm, respectively, for a field of view of 200 mm × 300 mm.


Particle position measurements

The single-shot methodology was first used to measure in real time the particle distribution and bubble dynamics in a cylindrical fluidized bed. Fluidized beds are vital to a variety of processes in the energy, chemical, and pharmaceutical industries and consist of a column filled with granular particles through which gas is introduced at the bottom. At sufficiently high velocities, that is, exceeding minimum fluidization velocity, particles become suspended, creating a fluid-like state. In many cases, bubbles of gas rise through the particulate phase in a fluidized bed. These bubbles promote mixing of particles and efficient heat transfer, making fluidized beds advantageous as compared to stationary beds of particles, yet they also allow gas to pass through the bed without intimate contact with particles, which can limit mass transport and chemical reactions. Thus, understanding the complex dynamics of bubbling is essential for the design and optimization of reactors. This understanding is currently limited because previous measurements of particle motion around bubbles have been restricted to 2D beds that are strongly influenced by wall effects (29).

Figure 2 shows images of particle distribution through a central vertical slice (10 mm thickness) with an in-plane resolution of 3 mm × 5 mm and a temporal resolution of 7 ms (see fig. S1 for pulse sequence details). The first row of Fig. 2B shows the complex dynamics of a smaller bubble accelerating first to coalesce with a larger bubble, followed by the combined bubble splitting up into two side-by-side bubbles. The second row focuses on the initiation of bubble coalescence and shows with unprecedented detail changes in shape that can be captured with 7-ms resolution. Figure 2C shows the detail of the changes in shape as two bubbles merge, which would be missed with lower temporal resolution. Figure 2D exemplifies the use of thresholding to distinguish between bubbles and the particle phase, and Fig. 2E shows that smaller bubbles are observed lower in the bed and near the walls, whereas larger bubbles are observed near the center and higher up in the bed, due to coalescence of rising bubbles. Movie S1 displays the complex bubbling dynamics at different gas velocities. Future studies can use this technique to decipher mechanisms behind complex phenomena, such as bubble splitting, and to explore the effectiveness of techniques, such as the insertion of internal objects, to optimize bubble characteristics.

Fig. 2 Real-time MRI of bubble dynamics inside a cylindrical 3D gas-solid fluidized bed.

(A) Schematic of the imaging setup. Particles in a cylindrical bed were fluidized by air injected through a perforated distributor plate at the bottom of the bed. The MR receiver array was placed concentrically around the fluidized bed. (B) Instantaneous MR images reveal details of bubble dynamics in gas-solid fluidized beds (see also movie S1). The coalescence of two bubbles and the subsequent splitting of the merged bubble following an indentation from the top surface of the bubble can be observed. (C) The high temporal resolution of 7 ms allows tracking the position (shaded areas) and detecting the velocities (lines) of the bubbles involved in the coalescence process. The speed of the trailing bubble increases significantly once it enters the wake of the leading bubble. (D) A threshold is applied to the MR images to determine the centers of mass and the equivalent bubble diameters of the bubbles. (E) Spatial distribution of bubble centers of mass (position of colored dots) and their equivalent circular bubble diameter De (color and size of the dots). Larger bubbles are found toward the central upper part of the bed. The scribbled area at the top of the bed demarcates the eruption region which was omitted from the bubble size analysis.

Particle velocity measurements

To understand granular dynamics, it is essential to measure not only particle distribution but also particle velocity. We combined the imaging shown in Fig. 1 with MR phase-contrast velocimetry to take snapshots of particle velocity in granular flow systems with an acquisition time of 22.5 ms, representing an increase in temporal resolution of four orders of magnitude as compared to previous measurements (58). Figure 3A shows the particle velocity field around an isolated rising bubble in a fluidized bed with corresponding streamlines (white), whereas Fig. 3B shows the dynamics predicted by potential flow theory (30).

Fig. 3 Measured and simulated granular flow around an isolated bubble in a fluidized bed.

(A) MR phase-contrast velocimetry of an isolated air bubble (white) rising in a fluidized granular system. The temporal and spatial resolution (Δx × Δy × Δz) of the acquisition was 22.5 ms and 3 mm × 7 mm × 10 mm, respectively. The colors denote the magnitude |vxy|, the black arrows denote the direction, and the white lines denote streamlines of the in-plane particle velocity, vxy. The highest particle velocities occur in the central wake region of the bubble (red area). (B) Numerical simulation of potential flow around a similarly shaped bubble assuming impenetrable bubble walls and irrotational flow. Considerable differences are observable around the top and along the sides of the bubble (see text for discussion). (C) Vorticity ωxy = ∇ × vxy of the measured flow field. Significant vorticity is observed in the wake region and side regions. A region of inverted vorticity can be observed around the top of the bubble, which is probably caused by particles “raining” down through the bubble. See also movie S2.

The flow fields show similar circulation patterns, but theory underestimates the particle velocity in the wake and overestimates the velocity around the upper exterior of the bubble. Davidson’s theory (30) assumes irrotational flow and an unchanging bubble shape; previous measurements of particle velocity around 2D bubbles (29) have shown significant vorticity in the wake but negligible vorticity elsewhere and thus supported the validity of this theory. Figure 3C shows the vorticity in the experimental flow field, calculated using ω = rot v. For this 3D rising bubble, significant vorticity is generated not only in the wake beneath the bubble but also along the entire exterior of the bubble. The significant differences between measurements and theory drive the need to develop theory that allows for rotational flow, changing bubble shape (see movie S2), and penetration of particles across the bubble boundary to capture this flow more accurately, which is vital to the understanding of mixing in fluidized beds.

Granular shock wave measurements

We can easily tune the MRI toolkit shown in Fig. 1 to study a wide range of problems in granular dynamics with open fundamental questions and important applications. Here, we show one example of tuning to create images with high signal intensity for regions of coherent rapid motion of particles and low signal for regions with incoherent motion (see Materials and Methods). We use this imaging technique to study the impact of an intruder on a granular packing and subsequent wave propagation of coherent particle motion (31), a phenomenon that has analogies to shock waves in complex fluid flow and is relevant to a variety of applications ranging from landing space probes to understanding seismic activity (32). Pioneering work in this field has used photoelastic techniques (33) which were limited largely to 2D systems. Recent work has also applied high-speed ultrasound imaging to measure shocks in dense suspensions (34) using interstitial water as the signal source for measurement, limiting largely this technique to wet granular systems. Here, we demonstrate the noninvasive measurement of the propagation of a granular shock wave in a dry 3D assembly of grains (Fig. 4 and movie S3).

Fig. 4 Response of a loosely packed granular material to the impact of a spherical intruder.

(A) Three selected frames of an MR image series (B) showing the impact of a heavy intruder into a granular packing (Δt = 53 ms). The leftmost image in (A) shows an area of significant signal increase following the impact of the intruder compared to preimpact conditions. This bright area propagates concentrically from the point of impact, indicating a wave of coherent particle motion and acceleration. The signal increase is followed by a signal attenuation [dark concentric area around the sphere in the central image of (A)] caused by incoherent particle motion. This incoherent motion relaxes to an immobile state in the course of approximately 0.5 s (see also movie S3). (C) 1D spatiotemporal plot (Δt = 2.3 ms) of the same event allows for a detailed analysis of the granular response dynamics. The yellow band corresponds to the white fringe in (B). The slope of this band determines the propagation speed of the spherical wavefront. Signal fluctuations (f = 30 Hz) throughout large regions of the granular material can be observed following the impact of the intruder, yet dampen out within 1 s. To provide an in-depth analysis of this phenomenon, a detailed parametric study is required, which is, however, beyond the scope of the current work.

The leftmost image in Fig. 4B shows a uniform signal through the granular bed before impact of the sphere because particles are stationary. The second image shows a white (high signal) partial circular region formed directly below the bed surface just after impact, indicative of a shock wave forming and propagating concentrically from the point of impact. The third image shows that after ~50 ms, the sphere has descended into the granular bed and the shock wave has propagated further into the bed. The fifth image shows that after ~150 ms, the shock wave has propagated through the entire bed and left a dark region around the impact point indicative of incoherent motion of particles in the region, which is reminiscent of the increase in entropy observed after the passage of a shock wave in a fluid (35). Figure 4C shows a 1D spatiotemporal plot of the signal resolved in the vertical direction (averaged in the horizontal direction, but still slice-selective). This plot shows a strong signal from the initial shock wave (yellow band), followed by subsequent smaller waves recurring at a near constant frequency and dampened in magnitude over time. In the future, this technique can be used to characterize how the velocity and magnitude of shock waves and subsequent dampened oscillations vary with different impact conditions and pressures confining the grains.


Here, we have introduced a new MRI methodology for studying granular dynamics, which increases acquisition speed by up to four orders of magnitude, enabling temporal resolutions comparable to optical measurements. To demonstrate the versatility of this method, we took snapshots of phenomena ranging from chaotic bubbling in fluidized beds to granular shock waves with contrast based on particle concentration, particle velocity, and the coherence of particle motion. The method has already created new insights into the significant vorticity surrounding bubbles in 3D beds, exposing the need to improve analytical theory in this area; researchers can also use measurements with this MRI methodology to test and improve computational models. For this work, we have engineered artificial particles with increased signal yield and lifetimes compared to agricultural seeds. However, in studies for which SNR and spatial resolution are less critical, the presented pulse sequences could also be used to image agricultural seeds. We are confident that future studies can apply this technique to address physical problems ranging from jamming (36) and segregation (37) in granular flows to mechanisms underlying glass transition (38) to the dynamics of suspensions (39).


MRI system

All measurements were performed on a 3-T whole-body human MRI system (Achieva 3.0T, Philips Healthcare) with a bore diameter of 600 mm. A gradient strength of 0.031 T m−1 and a gradient slew rate of 0.200 T m−1 s−1 were used in the imaging pulse sequences. A custom-built 16-channel RF receiver array was used for parallel signal detection.

Design and construction of the 16-channel RF detector array

Signal detection in parallel MRI is based on Faraday induction in multiple RF resonators that are tuned to the Larmor frequency of the probe nuclei in the static magnetic field B0 (40). The achievable parallel scan acceleration factor R and SNR of the array depend on the size and position of the receiver coils. In general, close fitting and conformal detector arrays provide the highest sensitivity and scan acceleration due to the high fill factor and the localized sensitivity profiles of the individual coils of the array (25). To determine the optimal size and positions of the receiver coils, finite element simulations (COMSOL Multiphysics and MATLAB) that solve the quasi-stationary Maxwell’s equations were performed. The array configuration was optimized with regard to coil size and position to yield maximal SNR homogeneity within the sample volume. The number of coils, 16, and the dimensions of the cylindrical sample volume (inner diameter, 210 mm; height, 300 mm) were kept constant in the simulations. In addition, the top and bottom end of the cylinder were excluded from coil coverage, allowing for unrestricted access to the system (for example, for the provision of fluidizing air). The final array consisted of 12 rectangular surface coils (129 mm × 105 mm), arranged in two rows around the outer surface of the cylinder, two circular coils at the top and the bottom of the cylinder, and two strip line coils (60 mm × 120 mm) at the cylinder edges perpendicular to the static magnetic field axis. Immediately adjacent resonators were overlapped geometrically by 10% of their area to reduce mutual inductance and resonance splitting (40). The array was realized on an acrylic tube onto which the 16 coils were glued using copper adhesive tape. Each coil was connected to a high-impedance preamplifier circuit that allowed for tuning each resonator to the Larmor frequency and matching the impedance of its port to the noise optimal impedance to maximize power transmission of the preamplifier (41, 42). To eliminate coupling of the receiver resonators to the transmitting resonator and to protect the sensitive preamplifier from high-power RF excitation pulses, the individual coils were Q-spoiled during high-power RF transmission using a PIN diode–switched LC tank circuit.


Spherical noncohesive core-shell particles with tailored mechanical and nuclear magnetic resonance (NMR) properties were designed and engineered with the goal of prolonging the effective transverse magnetization T2* and increasing the observable spin density of conventional granular signal sources (namely, agricultural seeds). The middle-chain triglyceride oil core of these particles makes up 73 weight % (wt %) of the total mass of a particle (see fig. S2). The core is contained in a 100-μm-thick shell of agar. Under experimental conditions (relative humidity of the fluidizing air was approximately 17%), the particles have an optically measured particle diameter of dp = 1.02 ± 0.12 mm, density ρp = 1040 kg m−3, dynamic angle of repose θr = 28 ± 2°, friction coefficient μ = tan(θr) = 0.54 ± 0.05, and coefficient of restitution e = 0.70 ± 0.03. Signal intensity and lifetime T2* were measured in a 3-T system by fitting an exponential function to the free induction decay. Compared to poppy and mustard seeds (T2* of 1.39 and 1.48 ms, respectively), the engineered particles yield a T2* of 1.75 ms due to their magnetostatically favorable spherical distribution of spins and an increased signal intensity due to the higher oil content. The knowledge of the NMR relaxation times T1 and T2* allows optimization of the flip angle θ = θErnst = cos−1 (e−(TR/T1)) (43) for any repetition time TR in order to maximize signal yield and to adjust the acquisition time Tacq = O(T2*) in order to determine the effective spatial resolution of the acquisition. The particles were fabricated by pumping the liquid precursor solutions of the core and shell through two concentric oscillating nozzles. Moisture can alter the mechanical properties of the particles through particle swelling and modification of their surface properties. Under experimental conditions, the particles contain 1.2 wt % humidity, which is stable over time periods of months if the humidity of the fluidizing air is stable and if the particles are stored in sealed bags between experiments. This moisture can be removed through vacuum drying. Only negligible electrostatic charging effects were observed in our experiments.

Fluidized bed

A clear acrylic tube (outer diameter, 200 mm; inner diameter, 190 mm; height, 300 mm) was flanged onto a perforated acrylic distributor plate. Below the distributor, a cylindrical wind box (outer diameter, 200 mm; inner diameter, 190 mm; height, 150 mm) was attached. Typical bed heights (unfluidized) were in the range of 170 to 240 mm, depending on the experiment. The top of the bed was open to the atmosphere. To ensure a homogeneous gas injection, the pressure drop across the distributor (10 mm thickness and 6416 laser cut holes of diameter 0.5 mm) was at least 10 times higher than the pressure drop across the fluidized bed at minimum fluidization. The homogeneity of the gas injection was improved further by injecting the air in the wind box via a perforated tube covered by 50 mm of foam rubber, a perforated plate, and another layer of foam rubber. The air flow was controlled with a mass flow controller (F-203AV, Bronkhorst).

Spin density particle position measurements

Instantaneous snapshots (acquisition time Tacq = 4.6 ms) of bubbling fluidized beds were acquired using a single-shot EPI readout with a temporal and spatial resolution of 7 ms and 3 mm × 5 mm × 10 mm, respectively (see fig. S1A for sequence details). Gas flow rates between zero and two times the minimum fluidization velocity umf (0.28 m s−1) were applied, and the unfluidized fill height of particles was 240 mm. Air bubbles (signal voids) were detected by applying a threshold of ρt = 0.43 times the mean particle signal intensity. The bubble areas Ab were detected and converted into their equivalent diameter of a circular bubble De=4Abπ.

Phase-contrast particle velocity measurements

The cylindrical bed described above was filled with particles to a height of 240 mm. A gas flow rate of umf was applied, whereas single air bubbles were injected into the bed through a nozzle (length, 40 mm; inner diameter, 6 mm) mounted at the center of the distributor plate. MR phase-contrast velocimetry was realized by applying bipolar motion-encoding gradients before the single-shot EPI readout (see fig. S1B for sequence details). These motion-encoding gradients yield a residual phase shift for particles with a nonzero velocity (26). To determine the two in-plane velocity components vx and vy, three data sets have to be acquired: one for each direction of motion and one reference data set without a bipolar velocity-encoding gradient. The temporal resolution of the acquisition is therefore roughly three times the repetition time, TR. The encoded field of flow was venc = 0.6 and 1 m s−1 for vx and vy, respectively. In our measurements (Fig. 3 and movie S2), we achieved a temporal resolution of 22.5 ms and an in-plane spatial resolution of 3 mm × 7 mm. The value of the spatial resolution in the phase-encoding direction, Δy = 7 mm, is coarser than the nominal value (5 mm), which assumes perfect Hermitian symmetry. This is due to the fact that the Hermitian symmetry is perturbed by sample motion that can lead to a reduction of the spatial resolution in the phase-encoding direction. Homodyne phase–constrained reconstruction (26) was used to reconstruct the velocity information for the acquired data (ReconFrame 3.0, GyroTools LLC).

Numerical simulations of particle velocity

Numerical simulations of granular flow around an air bubble were performed in COMSOL Multiphysics. The bubble shape was extracted from the MRI experiments (Fig. 3A). The model of Davidson (30) was applied, assuming (i) irrotational (potential) flow, (ii) incompressible flow, (iii) unchanging bubble shape, and (iv) zero penetration of particles through the bubble. Therefore, in this model, the velocity field can be described as the gradient of a velocity potential φ, which obeys the Laplace equation ∇2 φ = 0 globally. A cylindrical domain was modeled with free slip, zero penetration boundary conditions at its side walls and at the bubble surface. In the simulation, the bubble was kept at a fixed position in the domain, and instead, the granular material was set into motion by introducing a flux source that enters the domain from the top surface (∂φ/∂y= −vbubble). At the bottom surface, a vertical outflow was enforced by applying Dirichlet boundary conditions (φ = 0). Negligible effects were observed when the domain size (190 mm × 200 mm versus 380 mm × 400 mm) and mesh size (365,324 versus 136,106 tetrahedral elements) were varied. To compare the particle velocity field produced by the simulation with that from MRI measurements, the bubble rise velocity vbubble was added to all of the velocity vectors produced from the simulation.

Detection of granular shock waves

In these experiments, the sample container was filled with granular material up to a height of 170 mm. A random loose packing (packing fraction φ = 0.59) was restored before each impact experiment by fluidizing the granular material above umf and gradually defluidizing the bed at a rate of 0.003 umf s−1. The impacting object was a hollow plastic sphere of diameter 50 mm filled with MRI transparent zirconia particles (diameter, 200 to 300 μm), yielding a total intruder density of ρi = 3000 kg m−3. The intruder was dropped onto the center of the bed from a height of 50 mm. The impact of the intruder and the granular response of the bed were imaged in a vertical, central plane using a gradient echo pulse sequence (see fig. S1C). Shock and incoherent particle motion sensitivity were introduced by applying a short repetition time TR = 2.3 ms and echo time TE = 1.2 ms, applying a large flip angle θ = 3 × θErnst, and using a gradient and RF spoiling scheme, which does not balance the second-order gradient moments. Hence, the signal acquired in one readout is superposed by the signal from previous excitations. Therefore, the sequence is very susceptible to mechanical perturbations that alter the superposition pattern. If particles move coherently, the signal contributions superpose constructively, yielding an observed signal increase of +42% compared to stationary particles (white fringe in Fig. 4, A and B). If the particles move incoherently, the readout gradient acts as a diffusion-encoding gradient, which effectively dephases the spin coherence of the excitation, and thus, we observe diffusion-attenuated signal (−35%). Conceptually similar effects are exploited in MR diffusion imaging (44) and intravoxel incoherent motion imaging (45). Both 2D sequences with a temporal and spatial resolution of 52 ms and 2.5 mm × 2.5 mm, respectively, and 1D sequences with a temporal and spatial resolution of 2.3 ms and 2.5 mm, respectively, were acquired.


Supplementary material for this article is available at

fig. S1. Sequence diagrams and properties of the presented MRI pulse sequences.

fig. S2. Mechanical and NMR relaxation properties of the engineered particles.

movie S1. Particle position measurements.

movie S2. Particle velocity measurements.

movie S3. Granular shock wave measurements.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank S. Gross, M. Weiger, F. Hennel, B. Wilm, T. Schmid, S. Ivannikov, and C. von Deuster for scientific discussions and S. Wheeler for building the fluidized bed setup. Funding: This work was supported by the Swiss National Science Foundation under grant no. 200021_153290. Author contributions: A.P., K.P.P., and C.R.M. conceived the project. T.T. designed the artificial particles. A.P., D.O.B., and K.P.P. modeled and constructed the receiver coil array. A.P. performed the experiments and the analysis. A.P. and C.M.B. performed the numerical simulations. A.P., C.M.B., K.P.P., and C.R.M. wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Data are stored on a password-protected computer and backed up on a network attached storage drive located at the Department of Information Technology and Electrical Engineering, ETH Zürich. Data will be made available upon request.

Stay Connected to Science Advances

Navigate This Article