Probing complex geophysical geometries with chattering dust

The modern energy economy and environmental infrastructure rely on the flow of fluids through fractures in rock. Yet this flow cannot be imaged directly because rocks are opaque to most probes. Here we apply chattering dust, or chemically reactive grains of sucrose containing pockets of pressurized carbon dioxide, to study rock fractures. As a dust grain dissolves, the pockets burst and emit acoustic signals that are detected by distributed sets of external ultrasonic sensors that track the dust movement through fracture systems. The dust particles travel through locally varying fracture apertures with varying speeds and provide information about internal fracture geometry, flow paths and bottlenecks. Chattering dust particles have an advantage over chemical sensors because they do not need to be collected, and over passive tracers because the chattering dust delineates the transport path. The current laboratory work has potential to scale up to near-borehole applications in the field.


Supplementary
probability maps were exported as 2D tiff images, and then segmentation was performed on the full 3D stack of 2D tiff images.
The segmented images were then analyzed in DragonFly to calculate the volume, surface area, and center of mass (x,y,z positions) for each segmented bubble. Bubbles included in the distributions shown in Figures 1 & 2 in the manuscript and here in Supplementary Figure 1 were limited to those bubbles with volumes > 4πr 3 /3 where r = 5.0*pixel size. This requirement ensures diameters with at least 10 pixels. Figure 1 in the manuscript contains the results from high resolution imaging (20x objective) of a sub-volume from a single dust grain (sample G in Supplementary Table 1). From the 2D projection ( Figure 1a) and 3D reconstruction (Supplementary Figure 1b) a range of bubble sizes is observed. Using the interpretation condition (volumes > 4πr 3 /3 where r = 5.0*pixel size), only bubbles with radii greater than 4.5 µm were used in the analysis. For this limited field of view, individual bubble volumes ranged from 500 to 10 7 µm 3 , with the most probable bubble volume at this resolution being ~10 3 µm 3 (~ 6.2 µm radius for a spherical bubble).
For samples CDA1-CD6A, the entire grain was imaged with a 4x objective. Again, only bubbles with radii > 5*pixel resolution (assuming a sphere) were used in the determination of the bubble volume probabilistic distribution. This limits the interpretation of bubble volume to bubbles with radii greater than 21 to 26 µm. The spatial distribution of bubbles is shown in 3D plots of the center of mass for each bubble for samples CD2A & CD6A ( Figure 2 in the manuscript), CD3A & CD4A in Supplementary Figures 1c&f. For this field of view, both bimodal (e.g. Figure 2a in manuscript) and non-bimodal (e.g. Figure 2b) probability distributions of bubble volume were observed. For the 9 imaged grains, individual bubble volumes ranged between from 3.2 x 10 4 to 2 x 10 7 µm 3 , with the most probable bubble volume ~10 5 µm 3 (~ 28.8 µm radius for a spherical  Table 2 contains the number of bubbles and grain volume for each chattering dust grain. From these data, on average there are ~100 gas bubbles/mm 3 (Supplementary Figure 2).

Supplementary Note 2.
Methodology for relating acoustic hits to bubble size distributions. After 3-D X-ray microscopy, acoustic emission measurements were performed on the 9 individual chattering dust grains (Supplementary Table 1  After filling the cuvette with water, the AE sensors were connected to an AE measurement system (24 Channel Mistra Express) through preamplifiers (Mistra 1220-5054, 20/40/60 dB single-ended powered preamplifier) to record signals using Mistra AEWin software. The threshold amplitude for detection was set at 45 dB (with a 60 dB preamplifier setting with 100kHz -400Khz window) which was determined to eliminate ambient noise for this experimental set-up. An acoustic emission (referred to as a hit) was recorded when the signal amplitude exceeded this threshold. The amplitude of an AE signal is defined as the maximum positive or negative signal excursion during an AE hit. Supplementary Figure 3 shows selected signals emitted from sample grain CD04 during the course of the experiment from the 1 st recorded hit at time 17.8273442 seconds to a time of 302.9656898 seconds (Note: the AE system records experimental time out to the 7 th decimal place which is 0.1 µs). Each signal is normalized by its maximum positive amplitude to enable visualization of the waveforms. Peak amplitudes varied from 0.177 to 8.91 Volts (45 dB to 79 dB). The AE software converts Volts to dB: dB = 20 log (Vmax/1 µVolt) -(Preamplifier Gain in dB). The AEWin software determined the time and amplitude of every hit and saved the information to a summary line file. In addition, the signals were also recorded (2 MSPS, 100 µs pre-trigger). If more than one event is observed in a single waveform, only one hit per signal was counted (Supplementary Table 4). This potentially leads to an undercounting of AE hits.
Supplementary Figure 4 (b) shows the number of AE hits as a function of number of bubbles determined from the 3D X-ray microscopy (see Supplementary Note 1). The number of AE hits increases with an increase in the number of bubbles in a grain. A linear fit to the data shows that fewer events were recorded than the number of bubbles in a sample, which suggests that other events may have been below the selected threshold value for AE detection or occurred within the 100 millisecond window of another event and were not counted. Only for a bubble count > 3000, did the number of events come close to the bubble count.
The amplitude of the recorded hits on all of the channels ranged between 45 to 79 dB (with the lower limit set by the threshold).   The sphere with transducers was placed on a stand in a tank of water (300 mm x 300 mm x 425 mm). The immersion transducers were connected to the AE measurement system (24 Channel Mistra Express) through preamplifiers (Mistra 1220-5054, 20/40/60 dB single-ended powered preamplifier) to record signals (2 MSps, 200 µs pre-trigger, 2048 points, 0.5 µs/point) using Mistra AEWin software. The threshold amplitude for detection was set at 45 dB (with a 40 dB preamplifier setting with 100kHz -400Khz window) which was determined to eliminate ambient noise for these experimental conditions. The AE system was initiated and then a single dust grain was released through a hole at the top of the sphere and then fell under gravity. Supplementary Figure 5c compares signals from the 6 transducer pairs with each pair diametrically opposed and at the same depth, and the signal from transducer 1 at the bottom of the set-up. The 1 st motion at each sensor shows a positive excursion in amplitude indicating that the source is an explosive source.

Supplementary Note 4.
Methodology for relating dissolution rate and acoustic events. The bottom of a 100 mm diameter petri dish was attached with hot glue to a single AE sensor (Mistra F15A, F-series, Passive Wideband alpha sensor 100-450 kHz) to record the change in acoustic emissions as a particle dissolved in water. For these experiments, the sensor was connected to the previously described AE system with a preamplifier setting of 60 dB and a threshold of 25 dB which was sufficient to eliminate ambient noise for these experimental conditions. The AEWin software recorded the time and amplitude of each hit that was saved to a summary line file. The number of hits from each sample is listed in Supplementary Table 5. A custom-built digital imaging system was used to record images of a grain as it dissolved in water (e.g. Figure 3 in manuscript). The system consisted of a Spy camera (Spy camera for Raspberry PI No 1397 from Adafruit) connected to a Raspberry Pi Model B+ with 512MB RAM that captured 3 layered images (red, green and blue (rgb) arrays with dimensions of 2592 pixels x 1944 pixels) at a rate of 1 frames per second. The images were stored as png files directly onto a flash drive. The camera was mounted on a stand at a fixed distance from the sample to obtain images of the grain above the transducer with a resolution of a pixel edge length of ~55 µm. Image analysis was performed with a custom code written in IDL that performed image segmentation by taking the 1st layer of an image and dividing it by the 2nd layer. The location of pixels above a threshold value (see Supplementary Table 5) was used to identify the single dust grain (Figure 3a&d in manuscript blue regions). Pixel counting was used to find the area of a grain 2 . After initiating the AE monitoring and image capture systems, a single dust grain was placed in the water-filled petri dish with tweezers. The water was changed between tests. Of the 12 samples tested, only 3 samples (samples DT13, DT14, DT15 in Supplementary Table 5) disaggregated into 2 or more subparticles over time (e.g. DT14 in Figure 3d in manuscript). The functional form of the time rate of change of the cross-sectional area, ΔA/Δt, of the grains differed between particles that exhibited disaggregation and those that did not. Figure Supplementary Table 6 for transducers locations and see Figure 4b&c in manuscript) which were subsequently connected to the previously described an AE monitoring system (see section 2 but a system with only 8 channels and with a preamplifier setting of 40dB and a threshold of 40 dB which in this experiment was sufficient to eliminate ambient noise) to record the full waveform signal (2 MSPs, pretrigger 100 µs, 2048 points). Calibration experiments on location interpretation for the seeded fracture experiment were performed for the transducers location listed in Supplementary Table 6. A drop of water was placed on the rough surface of the fracture and then a single dust grain was released into the drop. Using differences in arrival time differences from among 3 or more transducers with amplitudes greater than or equal to 60 dB, the location interpreted from the AE signals was within + 5 mm of the location determined from an image of the dust grain on the fracture surface.
The 8 sensors were attached to the top half of the sample. The lower half was seeded with single grains of chattering dust. The upper half was then placed on the lower half, and a gentle pressure was applied to ensure the two fracture surfaces were registered. This process led to crushing of the dust, resulting in several dust fragments in some locations (locations with multiple stars in Figure 4 of the manuscript). Images of the initial seed location were taken for single surface and after registering with the lower half of the sample for comparison. The unsealed fracture sample was placed in a basin (42.9 cm x 35.6 cm x 18.7 cm or ~ 17 L) and inclined 6 o above the horizontal. The RaspberryPI image capture system described in section 4 of the Supplement Methods was set-up to record images of the fluid invasion during the course of an experiment. The camera was mounted at an angle a fixed distance from the fracture plane to capture the entire fracture plane, and a light was used to illuminate the fracture plane from below. The basin was white and acted as a light diffuser. Each image was recorded to a jpg file with an image size of 2592 pixels x 1944 pixels x 3 layers at a rate of 1 frame every 2 seconds. Using a constant head (3,731 Pa), water flowed into the basin and was then imbibed into the sample. The pixel edge length was ~160 micrometers. The image capture and AE systems were initiated once the water had imbibed into the lower right corner of the sample (as oriented in Figure 4b&c in the manuscript) and continued throughout the invasion processes (~ 25 minutes). Time syncing between the images and AE events is within +/-2 seconds. The first hit on the AE system occurred after 10 minutes which was the time for the water to invade to the location of the first dust grain.
Data analysis methods for results from seeded fracture experiments. Image analysis was performed using custom code written in IDL to determine the invasion front as a function of time. Supplementary Table 6. Positions of AE sensors for the coordinate system for seeded fracture experiment shown in Figure 4 in the manuscript. To identify the front, the difference between the 1st layer of an image and the 1 st layer of an image from 20 second range prior to the current image was taken. This provides an image of the regions invaded by the fluid in a 20 seconds. Using an image from earlier in the experiment (~56 seconds), a mask of the transducers is created from the 1 st layer of the image by finding all pixels less than a threshold of 80 which are then set equal to 1. A dilation algorithm was applied to the mask with a kernel to remove pixel jitter. Then the mask was inverted such that sensor locations were set equal to 0 while all other locations were set equal to 1. Each image was multiplied by the mask, after which a label_region command was used. The front was identified by finding the regions with the highest pixel counts and setting these regions to 1 and discarding the other regions. The color contours shown in Figure 4c in the manuscript represent a 20 second time window on the location of the front with 60 seconds between contours (edge-to-edge).

AE
Interpretation of the location of the chattering dust grain was based on the interpretation by the AEWin software and restricting the interpretation to events recorded by 4-8 sensors with amplitudes greater than 60 dB. To compare the time and location of an AE event with the fluid invasion front, the images were oriented (0.2 degrees counterclockwise rotation to an axis out of the plane) and scaled to match the AE coordinate system. The x-and y-dimensions of the images where scaled to account for the angle between the camera and the fracture plane (5 o and 19 o , respectively) and with a rigid shift of 5 mm in x and y (which is ~ 3.3% of the total length of the sample). The color of the symbols represents the time of the AE event which uses the same color scale as used to indicate the time of invasion for the contours in Figure 4b in the manuscript.

Uniform aperture and variable aperture fractures experimental set-up.
Transparent synthetic fractures were fabricated from two acrylic blocks measuring 150 mm x 150 mm x 100 mm (Supplementary Figure 8b) to enable video imaging of the descent of a grain under gravity in a fracture. The aperture was created by separating the blocks by a fixed distance. Uniform apertures of 0.5, 1, 2, 4, 8 and 10 mm were studied. For the variable-aperture fractures, 2 mm thick rubber was cut to create the in-plane variable geometry (Supplementary Figure 8c and Figure 6a-d in the manuscript) and was placed between the two acrylic blocks (Figure 6 in manuscript). 8 AE sensors (Mistra F15A, F-series, Passive Wideband alpha sensor 100-450 kHz) were attached to a sample with hot glue (see Supplementary Table 7 for transducers locations and example signals in Supplementary Figure 7) which were subsequently connected to the previously described AE monitoring system (see section 2 but with a preamplifier setting of 60dB and a threshold of 35 dB which in this experiment was sufficient to eliminate ambient noise) to record the full waveform signal (2 MSPs, pretrigger 100 µs, 2048 points).
The RaspberryPI image capture system described in section 4 of the Supplemental Methods was used to video the descent of each grain during each experiment. The camera was mounted a fix distance from the fracture plane to capture the entire plane. The video was recorded as a

Interpretation of chattering dust location from video images:
For analysis, the videos were converted using Free WMV AVI Converter software to 640 pixel x 480 pixels in with a frame rate of 25 fps, and then exported from QuickTime Player 7 as png images (640 pixel x 480 pixels). The pixel resolution in the images ~470 µm edge length. A custom code was written in IDL to segment the images to identify a dust grain and track the location of the dust as it fell. This was achieved by taking a subimage (50 pixels wide by 315 pixels in length) and then thresholding the image. Using the data cube (x-y-t, where t=time), two arrays were created from the threshold images to observe the lateral and vertical path from 20-500 images (# in sum depends on aperture and speed of descent): (1) a single image from the sum of y-t planes, and (2) a single image from the same data cube but of the x-t plane. These two arrays were then thinned, and the x-and ylocations of the maximum pixel value were found and exported to a text file along with the image time. The speed of descent was determined from a linear fit to the data in the graphing software Kaliedgraph.
Interpretation of chattering dust location from acoustic emission. Custom codes were written in IDL to (1)  A dust grain can float or remain at the air-water interface at the top of fracture for two reasons: (1) an initial hydrophobicity that maintains the particle at the air-water interface until it starts to dissolve; and (2) the particle dimension is too large to fit into the fracture aperture. The average float time as a function of aperture is shown in  (1). For apertures < 2mm, the float times are long because the particles were not small enough to fit into the aperture. Assuming that a spherical particle with a radius ~ 2 mm dissolves to 0.4 mm to fit within the 0.5 mm aperture, the time rate of change in cross-sectional area of a grain would need to be ΔA/Δt ~ -0.101 mm 2 /s (assuming a circular cross-section) for the average float time of 120 seconds (Supplementary Figure 9a). This value of ΔA/Δt is comparable to the value from the dissolution experiment described in section 4 of this supplement (ΔA/Δt -0.115 mm 2 /s with a standard error of 0.014 mm 2 /s). It is also important to note that dissolution is a function of surface area as dissolution rate decreases with decreasing surface area.
As a particle descends, or is transported through a system, its dimensions change because of dissolution. Descent times, t descent , ranged from ~15 to 1.5 seconds for fractures with uniform apertures of 0.5 mm and 10 mm, respectively (with average velocity of descent = <V descent > ~ 10 mm/s and 100 mm/s). Using the average ΔA/Δt = -0.115 mm 2 /s, the A of a grain would decrease by ~1.3 mm 2 for the 0.5 mm aperture and by 0.415 mm 2 for the 10 mm aperture fracture or roughly a change in radius, Δr, (assuming a circular cross section) of 0.741 mm and 0.234 mm, respectively. This change in the radius of a grain would affect V descent because the drag from the walls decreases with distance from the wall. An indirect measure of the effect of particle dimensions is shown in Supplementary Figure 9bc&d. Supplementary Figure 9c illustrates the effect of disaggregation on interpreting V descent and also shows the repeatability of V descent and suggests an effect of particle size, with the 3 rd particle falling the faster and having been subjected to dissolution the longest. Supplementary Figure 9d examines the effect of volume on V descent . The size of a chattering dust grain used in these experiments was measured against a ruler for three orthogonal directions. The V descent interpreted from the video images exhibits a first order dependence that shows a linear increase in velocity with size when the particle is much less than the aperture but decreases in speed as the particle approaches the size of the aperture. This is hypothesized to explain the large variation in velocity observed in the measurements for the 1 mm aperture (Supplementary Figure 9c). Knowledge of the dissolution rate and time enable estimates of changes in aperture. From the data for the 1 mm aperture fracture (specifically Test013), a velocity was selected by determining the average location of the first 100 points when a dust grain was floating and finding the velocity that yielded the minimum value (i.e. depth y ~ 75 mm in Supplementary Figure  8b). This approach yielded the full possible descent path from 0 to 150 mm depth. Once the velocity was determined, events were located by using a Broyden non-linear solver to fit the x & y positions of the source and the time to source for sets of transducer combination (in groups of 3). The error bars shown in Figure 5a in the manuscript are based on the averaged location from all of the groups for a given event. The error bars in Figure 5b in the manuscript are from the average of 7-10 tests performed for each aperture. The Broyden approach is related to Newton's method for finding function zeroes, but it has higher efficiency because it calculates the entire Jacobian up front rather than at each iteration. The signal from transducer 1 was chosen arbitrarily as the reference signal. The set of equations solved by the Broyden routine are based on the the following equation: where t s is the travel time from the source to transducer 1, the subscript t i represents transducer "i", v is the velocity, and the location is given by x, y, and z. The z location of the source was assumed to be at the center of the fracture plane, z s =0. The minimum in the value of the first 100 V= 2630 m/s was used to find the position of the source (x s , y s ) and t s for every event during the experiment.

Analysis of interpretation of chattering dust location with and without assuming refraction.
A numerical study was performed to determine the error in the interpretation of dust grain settling velocity that can occur from not considering refraction of the emitted acoustic signal as it is transmitted from the source through water and into the acrylic block. A custom IDL code was written to perform forward modeling based on Fermat's principle of path of least time for the refracted path (Supplementary Figure 8a). The coordinate system for the simulations is given Supplementary Figure 8b, and the sensor locations in Supplementary Table 7.
The fracture wall defines a planar interface between the water and the acrylic matrix. The plane of the wall (150 mm x 150 mm x-y plane in Supplementary Figure 8b) was discretized into 50 µm x 50 µm elements. For each transducer location, the time from the center of the transducer to each wall element, t tw , on the fracture wall was calculated using the speed of sound in acrylic (V acrylic = 2730 m/s) and a path length, r acrylic , of (se2) where h is the total height of the fractures (150 mm y-direction). The fracture height was discretized in 1 mm increments. The emission from a dust grain falling under gravity in water was treated as a point source. The travel time for a wave that traveled from the source to the fracture wall, t sw, was calculated using the speed of sound in water (V water = 1480 m/s) and the path length, r water , of (se3) where a is the distance of the source to the fracture wall. The travel path was taken as the path of least time between a source and a receiver, namely, when t tw + t sw was a minimum. The minimum distance occurs when sin θ 1 /sin θ 2 = V water /V acrylic . For comparison, the travel time along the direct path (red path in Supplementary Figure 10a) was also calculated assuming a velocity along the entire ray direct path of a fitted-V acrylic (based on locating the floating dust see Supplementary Note 6) which is the assumption that was used in the analysis of the experimental data. Supplementary  Figure 10 provides a comparison of the travel time from a descending source to receiver position 1 (sensor 1 in Supplementary Table 7) for water-filled fractures with an aperture 0.5 mm and 10 mm assuming both that the ray is, and is not, refracted. The difference in travel time between the two ray paths is shown in Supplementary Figure 10b for all simulated and measured apertures. For apertures b = 0.5 to 2 mm, the time difference is below the resolution of AE monitoring system (dashed line). The AE system records signals with 0.5 µs/point. A dt = 0.5 µs in acrylic yields a difference in path length of 1.365 mm. For b > 2 mm, the difference in arrival time increases.
The travel times between each source position and 8 sensor locations (4 on the z= -0.075 x-y plane and 4 on the z=+0.075 mm x-y plane with the locations used in the experiments, see Supplementary  Figure 8b) were calculated to determine the error in x location and interpretation of the speed of descent of a dust particle. The time differences (assuming sensor 1 as the reference as in the data analysis of the experiments) were used in the Broyden non-linear solver described above to determine the location of the source as it fell along the path (x,z)=(0,0) for y= -0.075 to +0.075 mm in 1 mm increments for apertures 0.5, 1, 2, 4, 8, and 10 mm.
As in the experiments, a fitted velocity, fitted-V acrylic , was determined that yielded a total travel path of the dust ~150 mm. For a grain falling along the (x,z)=(0,0), Supplementary Tables 8 lists the average and standard deviation in the x location of the dust and the ratio of the speed of descent based on including refraction, V dr , to the speed of descent assuming a direct path, V dd , between the source and receiver, as assumed in the experiments by fitting V acrylic = 2717 m/s for all apertures. The x location deviated from x=0 by 157 µm to 2.6 mm with increasing aperture and resulted in 1-3% under prediction of the speed of descent, V descent (Supplementary Table 8). When V acrylic is interpreted for each aperture (Supplementary Table 9), the under-prediction in V descent decreases and the interpreted x location is off from x=0 by 100 µm to 550 µm with increasing aperture.
In general, deviations in the interpretation of dust location and speed of descent increase with increasing aperture. For the experiments presented here, the error from assuming a direct path in the interpretation of speed of descent is on the order of 1-3% which is smaller than the error from experimental measurements from averaging the location from sensor groups for each event. Methodology for transport of chattering dust through inverted T fracture intersection. The intersecting fractures used in the flow experiments were fabricated from 3 acrylic blocks as shown in Figure 6 in the manuscript. The horizontal block (bottom) measured 150 mm x 150 mm x 50 mm in height. The left block measured 100 mm x 150 mm x 150 mm and the right block measured 50 mm x 150 mm x 150 mm. The surfaces of the blocks were smooth except for the bottom face of the right block which exhibited saw-cut roughness. The entire sample was submerged in a water tank (diameter ~ 0.88 m and height ~ 0.66 m). The aperture of the fractures were ~4-5 mm.

Supplementary
The top of the vertical fracture was unsealed to enable water to be pumped through the vertical fracture and also the horizontal fracture. The left & right sides (as shown in Figure 6 in the manuscript) each contained 3 ports that could be opened or closed individually during transport. Fluid was pulled or pushed through the horizontal fracture with a 265 Lph water pump attached to one side of the fracture and the other end of the pump open to the water tank.
15 immersion transducers (Olympus Immersion Transducers V303-SU with central frequency 1 MHz) were used to monitor the system. The transducers were connected to AE system through the preamplifiers (see section 2 in this supplement). The threshold amplitude for detection was set at 35 dB (with a 60 dB preamplifier setting with 100kHz -400Khz window) which was determined to eliminate ambient noise for these experimental conditions. The signals were recorded (10 MSps, with a 100 µs pre-trigger) and saved, along with a summary line file, using Mistra AEWin software. The locations of the transducers are given in Supplementary  Table 10 for the coordinate system shown in Figure 6 of the manuscript.
The AE system was initiated prior to the release of a single dust grain in the vertical fracture. The dust was released through a narrow tube (~ 3 mm diameter) into the vertical fracture by applying a stream of water from the nozzle of a water bottle. This ensured the descent of the particle into the fracture system. The location of release was varied during the experiment. Interpretation of the location of the chattering dust grain was based on the interpretation by the AEWin software and restricting the interpretation to events recorded by 8-15 transducers.