Ultrasound imaging velocimetry in particle-laden flows: counteracting attenuation with correlation averaging

Ultrasound imaging velocimetry (UIV) refers to the technique wherein ultrasound images are analysed with 2D cross-correlation techniques developed originally in the framework of particle image velocimetry. Applying UIV to opaque, particle-laden multiphase flows have long been considered to be an attractive prospect. In this study, we demonstrate how fundamental differences in acoustical imaging, as compared to optical imaging, manifest themselves in the 2D cross-correlation analysis. A chief point of departure from conventional particle image velocimetry is the strong variation in the intensity profile of the acoustic wavefield, primarily caused by the attenuation of ultrasonic waves in particle-laden flows. Attenuation necessitates using a larger ensemble of correlation planes to obtain satisfactory time-averaged velocity profiles. For a given combination of imaging and flow conditions, attenuation sets upper limits on volume fraction, penetration depth, as well as temporal resolutions that may be accessed confidently. This behaviour is demonstrated in two experimental datasets and is also supported by a modified cross-correlation theory. The modification is brought about by incorporating a lumped model of ultrasonic backscattering in suspensions into existing spatial cross-correlation analysis. The two experimental datasets correspond to two distinct particle-laden pipe flows: (1) a neutrally buoyant non-Brownian suspension in a laboratory-scale flow facility, wherein particle sizes are comparable to the ultrasonic wavelength and (2) a non-Newtonian slurry in an industrial-scale flow facility, wherein particle sizes are much smaller than the ultrasonic wavelength. We illustrate how and to what extent correlation averaging can counteract the adversity caused by attenuation. The work herein offers a template for one to evaluate the performance of UIV in particle-laden flows.


Introduction and scope
Particle image velocimetry (PIV) is a matured measurement technique and has become an integral part of the experimental fluid mechanics lexicon over the past three decades. At 1 3 56 Page 2 of 14 the heart of PIV algorithms lie spatial cross-correlation analyses, whose results are translated to velocity fields. Since the cross-correlation analysis can quantify motion in general, techniques developed for PIV have been used in several other contexts as well. One notable example is its application to ultrasound images, commonly referred to Ultrasound Imaging Velocimetry (UIV) or echo-PIV (Poelma 2017).
Unlike the technology developed for PIV, commercially available ultrasound imaging systems are meant for medical imaging, and not for fluid mechanics research. Thus, UIV finds its foremost applications in cardiovascular research. Nevertheless, UIV can be extended to study other flows as well. In particular, the potential to study multiphase flows is a well-acknowledged advantage of ultrasonics (Powell 2008;Takeda 2012;Poelma 2020). This work is geared towards the application of UIV to particle-laden flows.
Owing to differences in the instrumentation and working principles, there is a key difference between the "illumination" in UIV for particle-laden flows and conventional PIV. PIV is typically used to study single-phase flows seeded with tracer particles, and the laser pulses used therein are extremely powerful. Partially owing to this, the effect of illumination intensity variation across the field-of-view, along the direction of propagation, can often be safely neglected. Of course, this assumption may break down if excessive seeding is used or if a high concentration of fluorescent dye is additionally present in the flow (Crimaldi 2008). The assumption of uniform "illumination" intensity is also usually valid when UIV is applied to a single-phase flow which has been seeded with tracer particles. However, when applying UIV to particle-laden flows with moderate/high particle loading (where particles are the dispersed phase and not tracers), this assumption might not hold true any more, as acoustical energy can vary strongly with imaging depth due to attenuation (absorption and/or scattering of sound). Several studies have acknowledged the negative effect of ultrasound attenuation on velocimetry of suspensions, but few have discussed it in detail (Wiklund and Stading 2008;Kotzé et al. 2008;Ricci et al. 2012;Kotzé et al. 2016;Kupsch et al. 2019).
Here, we take a closer look at the performance of UIV for particle-laden flows, especially how attenuation appears in the processing chain. We do so by considering two distinct particle-laden pipe flows: one at laboratory scale with the flow of neutrally buoyant suspensions, wherein particle sizes are comparable to the ultrasonic wavelength; the other at industrial scale with the flow of a non-Newtonian slurry, wherein particles are far smaller than the ultrasonic wavelength. Experimental observations are complemented and supported by a modified cross-correlation theory for UIV in particle-laden flows. This modified theory is obtained by incorporating a lumped model of ultrasonic scattering in a suspension (Dash et al. 2021) in the cross-correlation theory developed for PIV by Keane and Adrian (1992). This allows us to form heuristic arguments on the performance of UIV for particle-laden flows.
This work is aimed at researchers who consider using UIV to study particle-laden flows. We intend to convince the reader on the existence of ceilings for particle loading, imaging depths, and temporal resolutions that can be accessed confidently using UIV. Using the modified cross-correlation theory and two experimental case studies, we provide a template on how one can assess UIV results in particle-laden flows. Specifically, we wish to show why and how ensemble correlation averaging is the best way to compute time-averaged velocity profiles. Being aware of the extents to which UIV returns acceptable results ensures that one can possess foresight on the accessible parameter space, prior to planning an extensive measurement campaign.
The manuscript is arranged as follows: In Sect. 2, we modify theoretical formulations for spatial cross-correlation analysis developed for PIV, by including characteristics of ultrasonic scattering induced by the presence of particles. After describing the two distinct experimental flow configurations in Sect. 3, we discuss the quality of the experimental results in Sect. 4. We end by summarizing our main findings in Sect. 5.

Extending cross-correlation theory to UIV for particle-laden flows
In this section, we incorporate a simple, lumped model of ultrasonic backscattering in a suspension, which we evaluated in an earlier study (Dash et al. 2021), into the cross-correlation theory developed for single-exposure, double-frame PIV systems (Keane and Adrian 1992). While Sect. 2.1 and Sect. 2.2 provide background information, Sect. 2.3 is novel and elucidates the performance of UIV in particle-laden flows.

Cross-correlation theory for single-exposure, double-frame PIV
Single-exposure, double-frame PIV is the most relevant mode from the perspective of UIV. While performing spatial cross-correlation of two single-exposure images, one is interested in the location of the peak representing the particle displacement. The height of this signal peak, R , should ideally stand out compared to random correlations, commonly referred to as noise peaks. Assuming a homogeneous distribution of tracer particles, no local variation of mean image intensities, as well as no local presence of strong velocity gradients, the ensemble average of the correlation of the intensity fluctuations, R , for a given fixed flow field, is given by: Equation (1) based on Keane and Adrian (1992) has several terms on the right hand side which contribute to the signal peak.
1. P is the number of correlation planes added. Conventionally, it is unity. This term is relevant for ensemble correlation averaging, which will be explained later. 2. I is the intensity of illumination. The exponent of two appears since two images are cross-correlated. 3. N I is the image density and is directly proportional to the number density of tracers. This number is typically optimized. Very low values would lead to poor PIV results while very high values might bring about speckle and/ or multiphase flow effects. 4. F I is the in-plane loss of correlation, and is inversely proportional to the in-plane displacement of particles. Values close to unity can be achieved with advanced interrogation techniques. 5. F O is the out-of-plane loss of correlation, and is inversely proportional to the out-of-plane displacement. These losses are caused by particles leaving the measurement volume and thus cannot be accounted for in the dataprocessing stage. 6. 2 00 F represents the particle image self-correlation and is related to the shape of the particle image. In short, the shape of the displacement-correlation plane will be the convolution between the displacement distribution and this term.
Time-averaged velocity results from PIV measurements can be obtained in several ways (see Meinhart et al. 2000, Figure 1 therein). Vector averaging (average of instantaneous velocity measurements) is commonly used in turbulent flows. An alternative is the correlation averaging technique wherein all the displacement-correlation planes are summed prior to obtaining the time-averaged velocity field. Correlation averaging was primarily introduced to improve velocity estimation in case of poorer signal-to-noise ratios (Santiago et al. 1998;Delnoij et al. 1999). Basically, by adding several correlation planes, the desired signal peak continues to add up and grow, while the undesired noise peaks cancel out. The presence of turbulence will result in broadening of the correlation peak (Kähler et al. 2006).

Interaction between ultrasonic waves and suspensions
We consider a lumped model which encapsulates the nature of ultrasonic backscatter in a suspension. While there are more rigorous models that express the interaction between acoustic waves and suspensions (Thorne and Hanes 2002), these (1) R = P I 2 N I F I F O 2 00 F have been primarily derived for dilute suspensions of sediment being probed by a single-element ultrasonic transducer. The lumped model allows one to evaluate arbitrary suspensions using linear array transducers (arbitrary beamforms). This lumped model consists of two terms: the peak backscatter amplitude, I 0 , and the amplitude attenuation rate, I . The first term represents the acoustic energy backscattered by the initial layer of particles of a uniform suspension. The second term represents the decay of acoustic energy per unit length in the axial direction, for waves propagating through a uniform suspension. The magnitude of these terms typically increases with increasing volume fraction, , of the suspension. For simplicity, we disregard multiple scattering (Tourin et al. 2000;Snieder and Page 2007). This phenomenon nullifies the assumption that the time of a received signal can be converted to a corresponding depth. Moreover, we do not include propagation in the acoustic coupling medium, the wall materials, and specular reflections at walls (Wang et al. 2003). Including the role of other factors such as measurement volume, image pre-processing (envelope detection, log-compression), electronic time gain compensation would require a detailed, quantitative description of the fundamentals of acoustics and ultrasound imaging electronics.
Assuming a general, heterogeneous distribution of particles, (z) , across the imaging depth axis, z, the generalized backscatter profile can be defined by Eq. (2). The ultrasound transducer is located at z = 0 . This is analogous to equations 7, 8 in Thorne and Hanes (2002) and equation 5 in Saint-Michel et al. (2017). See Fig. 1(a) for the corresponding coordinate system. This means that the intensity of the ultrasonic wave corresponding to the depth z is composed of two parts. The prefactor is directly related to the volume fraction of particles at depth z. The exponential term, on the other hand, is related to the distribution of particles between the transducer and depth z.
The above equation resembles a Beer-Lambert law. Typical Beer-Lambert laws are based solely on attenuation of beams between a separate transmitter and receiver where the prefactor is the initial energy deposited (constant). On the other hand, ultrasound imaging (pulse-echo) is based on backscatter, with the same instrument functioning as a transmitter and receiver. Thus, the prefactor is a depth-dependent variable.

Cross-correlation theory for UIV of particle-laden flows
Due to the many differences between UIV for particle-laden flows and conventional PIV, several changes can arise in Eq. (1) by incorporating Eq. (2).
(2). 2. N I : This term is now dictated by the particle-laden flow itself. In particle-laden flows, the dispersed phase is generally used for velocimetry and no additional tracer particles are introduced. 3. 2 00 F : Ultrasound images of particles are sensitive to several factors. For example, the image of a point target is sensitive to its location in the insonified region (see Ortiz et al. 2012, Figure 2 therein). Particle image shapes are also determined by the relative wavelength defined as 2 a∕ , where a is the particle radius and is the ultrasonic wavelength (see Baddour et al., 2005, Figure 4 therein). 4. P, F I , F O : There is no major change.
For the following discussion only, we assume to be invariant in z such that particles are uniformly distributed with a bulk volume fraction b . Moreover, we introduce the exponent for ensemble size P, whose significance will be shown later. Typically, is unity. Combining all of the above yields Eq. (3). (3) As a next step, we take a logarithm on both sides of the above equation, and assume F I and F O to be unity, aiding the physical interpretation. This yields: The following discussion on how different factors contribute to the displacement-correlation peak ( R ) serves as a starting point for explaining the quality of the experimental data.
1. P ↑ ⟹ R ↑ : By summing more correlation planes (term A), the height of the peak should rise. 2. z ↑ ⟹ R ↓ : By going deeper into the imaging plane, more acoustic energy is lost due to attenuation. As a result, the signal-to-noise ratio of the displacementcorrelation plane may drop (term C). In flows with low image density (term D), the particle self-correlation term (term E) may also be dependent on z. 3. b ↑ ⟹ R ↕ : The effect of adding particles into the flow is ambiguous. Increasing b usually increases the magnitudes of I 0 and I , which have opposing effects Fig. 1 a Coordinate system for ultrasound imaging and global overview of experimental setups; b laboratory scale; c industrial scale on R (terms B, C). Since term C also includes z, increasing b is likely to create a stronger gradient in the quality of the cross-correlation planes along the imaging depth. Increasing b also increases N I (term D) which contributes positively to R . Thus, for a fixed a and , there shall be a critical b wherein the imaging will switch from particle-image to speckle mode (term E). 4. a (for fixed , either fix b or N I ): The effect of modifying the particle size on R is difficult to generalize. Increasing a increases the relative wavelength, 2 a∕ , which can have a strong effect on I 0 , I , and the particle image shape (terms B, C, E). The magnitude of both I 0 and I are expected to rise with increasing a, especially when N I is held constant. When a is increased while holding b constant, N I reduces (term D) in addition to effects on other terms as well (terms B, C, E).
In such a scenario, there will be a critical a wherein the imaging will switch from speckle to particle-image mode. 5. : The ultrasonic wavelength primarily influences the displacement-correlation plane via I 0 , I , and the particle image shape (terms B, C, E). Since the imaging resolution is dependent on , it could also affect N I (term D) and thus the imaging (particle-image/speckle mode). Generally, a larger wavelength leads to better penetration (reduced attenuation). Scharnowski and Kähler (2016) investigated the effect of spatially-invariant image noise on PIV results. They do so by formulating a general approach to quantify the effect of image noise from autocorrelation planes. To simulate the effect of image noise in experiments, the authors reduce the power contained in the laser pulses. In comparison to the above work, there are a couple of differences in our approach. For starters, the illumination intensity in UIV is non-uniform and depth-variant. Moreover, we directly incorporate the effect of noise into the intensity illumination term already present in Eq, (1), instead of using the autocorrelation.

Experimental datasets: particle-laden pipe flows
With the modified cross-correlation theory for UIV in particle-laden flows being established in Sect. 2, we now consider two experiments. One of the experiments is at a laboratory scale (total system volume ∼ 10 litres, pipe length ∼ 3 meters), while the other is at an industrial scale (total system volume ∼4000 litres, pipe length ∼ 85 meters). Shown in Fig. 1 are schematics of the experimental facilities, while key differences are presented in Table 1. The two experimental datasets are used to demonstrate that the modified cross-correlation theory sufficiently explains the quality of the experimental data, despite the differences in the experiments. While these two examples represent two distinct scales of experimental setups, the examples also represent two distinct scales from the perspective of interaction between ultrasonic waves and suspensions. The relative wavelength ( 2 a∕ ) in acoustics is analogous to that in optics with three regimes: long-wavelength / Rayleigh scattering ( 2 a∕ ≪ 1 ), intermediate-wavelength/Mie scattering ( 2 a∕ ≈ 1 ), shortwavelength/geometrical scattering ( 2 a∕ ≫ 1).

UIV at laboratory scale
The laboratory-scale facility was constructed to investigate the laminar-turbulent transition in neutrally buoyant non-Brownian suspensions. A precision glass pipe with an internal diameter of 10.00 mm is laden with suspensions of varying bulk volume fractions. The carrier liquid is a saline water mixture and the dispersed solids are rigid, spherical polystyrene spheres with diameters of 530 ± 75 m. An ultrasound imaging system (Verasonics Vantage 128) collects raw data for UIV at high temporal resolution (400-4000 frames per second). The system is coupled to a linear array transducer (L11-5v) of 128 elements with a pitch of 0.3 mm and an elevation focus at 18 mm. Ultrasound frequencies of 10.5 MHz are used. Data are sampled at 62.5 MHz with a 14-bit resolution. Image acquisition is performed in parallel, meaning all elements transmit ultrasonic pulses and later receive echoes simultaneously. Four sets of 10000 images are collected. At the location where ultrasound images are recorded (270 pipe diameters downstream of the pipe inlet), the glass pipe is replaced by a plexiglass pipe (internal diameter = 10 mm) to reduce reflection losses caused by the higher acoustic impedance difference between glass and water. Acoustic coupling between the transducer and the pipe is obtained via a water filled container surrounding the pipe. The volumetric flow rate is determined by collecting a sample of the suspension in a beaker (total volume collected per unit time).

UIV at industrial scale
The industrial-scale facility (" -loop", Deltares) was constructed to investigate the hydraulics of non-Newtonian slurries. A PVC pipe with an internal diameter of 100 mm is laden with a non-Newtonian slurry. Kaolin(ite), a form of clay, is mixed in tap water to form the slurry (21% w/w or 9% v/v). The particles are fine, with a median size of 5.18 m. Raw data for UIV were acquired using sequential beamforming techniques with a SonixTOUCH Research (Ultrasonix/ BK Ultrasound). This imaging system is coupled to a linear array transducer (L14-5/38) of 128 elements with a pitch of 0.3 mm and an elevation focus at 16 mm. Acoustic coupling between the ultrasound sensor and the pipe is achieved with Aquasonic 100 gel (Parker laboratories, Fairfield, NJ, USA). A custom "aligner" is 3D-printed, which allows for an accurate placement and alignment of the ultrasound transducer with respect to the pipe wall. The transducer is placed in the backward sweep mode (beam sweep direction opposite to flow direction) to obtain reduced particle displacements (Zhou et al. 2013

Data processing
For both sets of experiments, B-mode images are generated from the raw data. For the industrial-scale experiments, the time-averaged RF image is subtracted prior to envelope detection and log-compression to generate B-mode images.
In the laboratory-scale experiments, the B-mode images are first generated following which the time-averaged B-mode image is subtracted. For the laboratory-scale experiments, image sizes are 128 × 200 pixels (lateral × axial direction), whereas for the industrial-scale experiments typical sizes are 33 × 812 pixels. The B-mode images are processed using a slightly modified version of PIVware (Matlab script based on Westerweel 1993). We perform single-pass processing without any advanced techniques such as window offsets (Westerweel et al. 1997) or iterative image deformation (Scarano 2001). This simplicity is maintained so that results can be directly compared, without influence of these steps. For the laboratory-scale experiments, a grid of interrogation windows with sizes of 32 × 8 pixels are used with 50% overlap in both directions, resulting in six columns of velocity vectors. This grid is generated only between the pipe walls. For the industrial-scale experiments, window sizes of 16 × 24 pixels are used, resulting in three columns of velocity vectors. An overlap of 50% is used in the axial direction and 85% in the lateral direction. The 85% overlap is a practical necessity, given that the image size is 33 pixels in the lateral direction while the window size is 24 pixels. Doing so ensures three columns of vectors, one of which is not at the edge. The displacement-correlation peak is estimated with subpixel accuracy by fitting a 2D Gaussian function in a 5 × 5 pixel neighbourhood of the maximum. This choice is warranted due to peak broadening, as shall be shown later. Moreover, for the industrial-scale experiments, sweep correction is warranted (Zhou et al. 2013). While UIV is a 2D2C technique, for the purpose of this paper, we focus solely on the variation of the streamwise velocity along the wall-normal direction (axial direction of ultrasound images). In both cases, we pick a single column of vectors for further analysis.
Examples of instantaneous B-mode images from both experiments are shown in Fig. 2. To obtain the mean intensity profile, the total ensemble of images (threedimensional array-lateral direction, axial direction, time) is averaged along time and lateral direction. Even though the volume fraction of the dispersed phase is 0.09 in both examples ( Fig. 2(a,b)), the source density in (b) is several orders of magnitude higher than in (a), owing to the finer particle sizes. As a heuristic, the source density may be approximated by the concentration of scatterers per mm 3 (Poelma 2017). In both examples, attenuation of the ultrasound signal from top to bottom is evident. Due to the differences in the relative wavelengths, it can be expected that attenuation of the ultrasonic waves is Another visible difference is the imaging mode itself. In the laboratory-scale experiments, features of individual particles are discernible at smaller imaging depths, before the image devolves into speckle-like patterns at larger imaging depths due to a mix of attenuation and multiple scattering. In contrast, for the industrial-scale experiments, the slurry is visible as speckle, due to high source density.
For the laboratory-scale experiments, individual particle images are rather peculiar and the spatial extent of a single particle image far exceeds its physical dimensions (Fig. 2(c)). The shape of the particle images has clear artefacts. The topmost reflection likely corresponds to the top of the particle, whereas the others are "copies" of the same particle. For the single particle highlighted in Fig. 2(b), the spatial dimensions of the particle image are approximately 1.8 mm × 2.9 mm. Similar observations have been made in the past, and the interested reader is directed to Section V.A. of Baddour et al. (2005), wherein this observation is discussed in depth. One of the explanations offered for such behaviour is the scattering regime (high relative wavelength). A negative consequence of such particle images is that it may lead to inaccuracies in the UIV analysis (via term E in Eq. (4)).

Quality of cross-correlation results
Next, we view the quality of the cross-correlation results in the experiments through the lens of Eq. (4).

Effect of P, z with fixed b
We first consider a single example from the industrial-scale experiments wherein the loop is operated at a bulk velocity of U b = 1.48 m/s. This corresponds to the turbulent flow regime and the slurry is homogeneous, i.e. is not a function of z. Moreover, the experiments are performed for a single b . It must be noted that this velocity is higher than what is typically studied using UIV ( ∼ 0.5-1.0 m/s).
Corresponding correlation planes as a function of P and z are shown in Fig. 3. Expectations on the basis of Eq. (4) are qualitatively visible here. For example, the clarity of the displacement-correlation peak increases with increasing ensemble sizes and decreasing imaging depths. At the same time, this hints at the need for tens of thousands of images in order to quantify the time-averaged velocity at the pipe axis (around imaging depths of 55 mm). Moreover, with increasing imaging depths, the influence of noise grows, and can impact the accurate subpixel estimation of the displacement peak. Furthermore, if it is desired to extract Reynolds stresses from ensemble correlation averaging of UIV data (Scharnowski et al. 2012), the effect of noise would also need to be considered.

Fig. 3 Evolution of correlation maps for the industrial-scale experiments as function of P and z
For the purpose of quantifying the quality of the correlation result, we use the metric "peak-to-root-mean-square ratio". We refer to this as "signal-to-noise ratio (SNR)" henceforth. While PIV traditionally relies on the ratio between the two tallest peaks, a change is warranted for the present study. For starters, the peak ratio is considered an ad-hoc measure (Xue et al. 2014), and is best applicable for cross-correlation of two images (as in vector averaging). Moreover, since the quality of ultrasound images is inherently low, typical peak ratio values are < 2. Since we add several planes, SNR, as defined above, provides a better insight on its quality.
In Fig. 4(a), we plot the SNR of the correlation planes as a function of P and z. This serves as a quantitative counterpart to Fig. 3. Thus, it confirms that the quality of the cross-correlation plane improves with increasing ensemble size and decreasing imaging depths. The SNR seems to be best around imaging depths of 20 mm, while one would expect it to be better at small imaging depths. There could be a few reasons for this. The region close to the transducer falls in the "near-field" of the transducer (transmit focus is at 45 mm), where pressure amplitude characteristics of the ultrasonic waves are complex and unfavourable (see Bushberg and Boone 2011, Figure 14-14 therein). From the perspective of the flow, the presence of velocity gradients could also be detrimental to the SNR at imaging depths close to the pipe wall.
This plot, however, also allows one to judge whether time-resolved measurements are possible using sliding averaging (discussed later). For example, at an imaging depth of 20 mm, about 10 correlation planes are needed to obtain a SNR of approximately 6. In comparison, at an imaging depth of 40 mm, nearly 200 correlation planes are needed to obtain a SNR of approximately 6. This implies that in order to obtain identical quality of correlation planes, about 20 times more images would be required at 40 mm than at 20 mm. This is an important factor to consider prior to planning time-resolved measurements in unsteady flows with sliding averaging.
Previously, in Eq.
(3), we introduced the exponent for the ensemble size, P. For a stochastically stationary flow and with low levels of image noise, the exponent should be unity (see Adrian and Westerweel 2011, Figure 8.27 therein). This would mean that adding ten correlation planes would increase the peak height tenfold. In Fig. 4(b), we investigate whether this indeed holds true for experimental data. The contour plot therein shows the logarithm of the maximum value of R . As expected, this value increases with increasing P (Eq. (4)). To estimate at each imaging depth, a linear fit is performed between log 10 (max(R )) and log 10 P . The slope of the linear fit is .
The exponent is the growth rate of the maximum value of the correlation plane, which is also plotted. For imaging depths up to approximately 40 mm, is nearly unity, as expected. However, its value then gradually drops to a value of 0.5 with increasing imaging depths, after which it stays nearly constant. We speculate that the value of 0.5 is associated with additive properties of random noise. This is supported by Figure 8.27 in Adrian and Westerweel (2011), wherein the noise peak grows with an exponent of 0.65 (a value lower than unity). The depth at which the exponent becomes constant coincides with the depth where no velocimetry results can be derived. This critical penetration depth is indicated by a horizontal black line in Fig. 4. The short region where the exponent is 0.5 < < 1 might represent regions where noise starts playing an increasing role. For example, the received ultrasonic backscatter corresponding to these imaging depths could be pushing close to the resolving power of the electronics, subsequently leading to lower "pixel fill ratios" and quantization errors. As a result, Fig. 4 a SNR of summed correlation planes as function of P and z. b Maximum peak height as function of P and z. The panel on the right shows as a function of z. Horizontal black lines represent the critical penetration depth for the maximum number of images analysed increasingly lower fractions of correlation planes would contribute constructively. Thus, the situation could be salvaged by simply recording even more images. As mentioned previously in Sect. 2.1, there are numerous ways to quantify the time-averaged velocity profile. Here, we consider three different ways: vector averaging, sliding correlation averaging, and ensemble correlation averaging. Sliding correlation averaging represents a hybrid of the other two techniques, which allows attaining reasonably good signal-to-noise ratio while having a certain temporal resolution (Samarage et al. 2012). In this example, a total of 12514 cross-correlation planes are available. The ensemble correlation average is thus based on the sum of all planes. Vector averaging is based on taking the average of instantaneous velocities evaluated by interrogating images 1-2, 2-3, 3-4 ⋯ . Sliding correlation averaging is based on taking the average of "instantaneous" velocities evaluated by summing correlation planes 1-10, 6-15, 11-20 ⋯ . This selection of summing 10 planes is based on the trade-off between obtaining acceptable temporal resolution and acceptable signalto-noise ratios. Please note that 10 is not the optimal value but a typical one.
We compare the results of time-averaging from the above three techniques in Fig. 5(a). We expect that correlation averaging yields the most accurate results on basis of Fig. 3 (since the peak on the correlation plane begins to stand out as more planes are used). Evaluating the integral ∫ 1 0 2(u(r∕R)∕U b )(r∕R)d(r∕R) returns the value of 0.964 (expected value = 1), ascertaining the good performance of the correlation averaging technique. It can be seen that vector averaging as well as sliding averaging under-predict the velocity. Vector averaging, especially, performs poorly and is unreliable. Sliding correlation averaging returns lot more satisfactory results until about 33 mm deep, beyond which it fails. Correlation averaging works well until imaging depths of approximately 50 mm before failing. Of course, collecting even more images might push the penetration depth a bit further. Results for sliding averaging in Fig. 5(a) have been improved by neglecting detected outliers while computing the temporal mean. Further explanation is provided in Appendix A.
Thus, one must resort to ensemble correlation averaging rather than sliding or vector averaging to compute timeaveraged velocity profiles in particle-laden flows. The more correlation planes used, the deeper one can obtain velocities. Sliding averaging can retrieve time-resolved velocities, albeit until smaller imaging depths.

Effect of b
In order to investigate the effect of particle loading, we consider data corresponding to the laboratory-scale experiments, where b was varied from 0 to 0.20. For b ≥ 0.175 , we were unable to extract any kind of useful velocity information, effectively restricting our analysis to 0 ≤ b ≤ 0.15.
In Fig. 5(b), we plot time-averaged velocity profiles for a laboratory-scale experiment in a manner similar to Fig. 5(a). The profiles are also based on 12515 images. Sliding averaging already fails at imaging depths of approximately 3 mm, implying that any time-resolved measurements will be restricted to the near-wall region. Correlation averaging is able to salvage time-averaged velocities until depths of 9 mm. Even though the volume fractions in Fig. 5(b) and Fig. 5(a) are similar, the limiting penetration depths are strikingly different.
Several factors contribute to the vastly differing penetration depths (based on Table 1). For example, the degree of specular reflection will be different because of different materials of the pipe wall, and will be dependent on the acoustic impedances. Differing material and thickness of  Fig. 7(d) the pipe wall would lead to different degrees of acoustic energy losses by absorption. Similarly, when using a parallel acquisition system, as compared to a sequential one, relatively inferior quality images are obtained in exchange for higher frame rates. Lastly, the different relative wavelengths result in different attenuation rates. Based on Figure 1(b)-(c) in Thorne and Hurther (2014), it can be expected that I 0 and I are higher for higher relative wavelengths. Thus, it can be expected that if all other factors remain unchanged, performing UIV on particle-laden flows with finer particles would return better results.
In Fig. 6, we plot the signal-to-noise ratio of the correlation planes as well as for a range of particle volume fractions. Figure 6(d) is a low velocity measurement (laminar conditions) wherein radial migration of particles is expected, leading to a non-uniform distribution of particles (see Dash et al. 2021, Figure 13 therein). The remaining cases have higher velocities and it can be assumed that particles are relatively uniformly dispersed. To facilitate a fair comparison, an attempt was made to hold F I constant, which is directly related to the particle displacements. For example, Fig. 6(d) corresponds to measurements in a flow with U b = 0.06 m/s and images were recorded at 1000 frames per second, while for Fig. 6(e), U b = 0.59 m/s and images were recorded at 4000 frames per second. Thus, for the former, image 1 is correlated with image 10, image 2 with image 11, image 3 with image 12 and so on. For the latter, image 1 is correlated with image 5, image 2 with image 6, image 3 with image 7 and so on. This gives comparable pixel displacements across the cases.
The effect of varying bulk volume fractions is visible in Fig. 6(a-c, e-f). For b = 0.02 , the SNR is nearly uniform across the entire P − z plane, suggesting that time-resolved measurements with depth-invariant quality can be taken, implying that attenuation, Iz , does not play a role yet (term C, Eq. (4)). When b is increased to 0.04, the SNR improves uniformly due to increase in I 0 (term B, Eq. (4)). Upon increasing b to 0.06, the SNR improves slightly, at least for z < 8 mm due to I 0 , but beyond this, there is a drop in SNR due to the growing dominance of Iz . For b = 0.08, 0.15 , the conflicting roles of I 0 and Iz cause a starker contrast between smaller and larger imaging depths. Only for b = 0.15 , the exponent shows a gradual decline from 1 to 0.5, like Fig. 4(b).
For simplifying the interpretation of Eq. (4), it was assumed that the particles were uniformly distributed throughout the flow. An example where this might not hold true is when particle migration occurs. At low Reynolds numbers for semi-dilute suspensions, inertial migration has been reported (Ardekani et al. 2018;Dash et al. 2021;Yousefi et al. 2021). We compare SNR patterns in Fig. 6(d-e) with respect to the presence/absence of radial migration for constant b . When radial migration occurs, there is locally higher particle volume fraction between the pipe axis and wall, leading to high I 0 and increased SNR. However, this is also accompanied by increasing influence of I , leading to a quicker decay in SNR values. Interestingly, Fig. 6 SNR of cross-correlation planes and as function of P and z for laboratory-scale experiments the signature of inertial migration, with reduced volume fraction at the pipe axis, is weakly reflected in the SNR contour plots as well.

Conclusions and Outlook
We elaborated upon the differences that arise between UIV for particle-laden flows and conventional PIV. Cross-correlation theory developed for PIV is modified to include ultrasonic backscatter behaviour in suspensions. This allows us to heuristically argue how different factors propagate into the cross-correlation plane of UIV measurements. This modified theory sufficiently explains the experimental data and can be seen as a step towards a complete understanding of the role of ultrasound imaging in UIV for particle-laden flows.
The key takeaway message is that while applying UIV to particle-laden flows, one must be aware of the presence of a ceiling for particle volume fraction, imaging depth, as well as achievable (depth-variant) temporal resolutions. These ceilings are mainly imposed due to the attenuation of ultrasound waves. For obtaining time-averaged velocities, ensemble correlation averaging is strongly advised. It is strongly recommended to work at the level of correlation planes while assessing the quality of velocimetry results. Determining (growth rate of correlation peak height as function of number of images) on a preliminary measurement can help one assess whether recording even more ultrasound images will be beneficial or not. Sliding averaging can help to gain insight in time-varying flows. Naturally, the depth at which it will become unreliable is reached sooner than using correlation averaging of the entire data set. From a practical perspective, combining UIV with finer particles is likely to be more rewarding than with coarser particles.
Our study highlights the need for one to perform preliminary tests in an ex-situ mock-up prior to an exhaustive and resource intensive measurement campaign. Even though there is no harm in performing measurements and realizing the limitations a posteriori, valuable time and labour can be saved while possessing realistic expectations of the feasibility a priori. This is especially valid when large volumes of suspensions need to be prepared. If the generation of realistic, synthetic ultrasound images for particle-laden flows can be realised (similar to Monte Carlo simulations for PIV), it will greatly help assess the performance of UIV in a controlled manner. As an example, Swillens et al. (2010a), Swillens et al. (2010b) perform speckle tracking on synthetic ultrasound images of cardiovascular flows.
In conclusion, we provide a methodology and framework (for example, Eq. (4), SNR, ) for one to evaluate particleladen flows with UIV. Despite the few limitations of UIV in the context of particle-laden flows, it enables acquiring unique and useful experimental data.

A Quality of sliding averaging
Previously, in Fig. 5(a), sliding correlation averaging was used to compute a time-averaged velocity profile. Here, we present how the quality of this time-averaged velocity profile was improved. Shown in Fig. 7(a)-(c) is an example wherein the detection of outliers for sliding averaging technique is demonstrated. While typical outlier detection / replacement techniques take place in spatial coordinates (Westerweel and Scarano 2005;Garcia 2011), in Fig. 5(a) we have performed this analysis on spatio-temporal coordinates. Such a change is necessary since only three columns of velocity vectors are computed, with only the middle one being reliable, ruling Fig. 7 Space-time plots of a streamwise velocity, b radial velocity, c SNR corresponding to Fig. 5(a). The white dots indicate points where SNR < 3 . d Joint probability distribution function between signal-to-noise ratios of displacement-correlation planes and "instantaneous" streamwise velocity measurements at imaging depth of approximately 12 mm (horizontal dotted line in Fig. 5(a)) out post-processing in spatial coordinates. In short, for each point on the spatio-temporal plot, the SNR of the corresponding correlation plane is evaluated. If this SNR < 3, we do not include it while computing the mean. For imaging depths > 35 mm, it can be seen that the SNR based on 10 planes is simply unacceptable. We do not attempt to replace the detected outliers. Replacement strategies will function well at low velocities, as long as outliers are not clustered together. At higher velocities, the higher fraction of outliers coupled with clustering of outliers will overwhelm any correction procedure. For example, at 19.9 s between depths of 30-40 mm, the few detected outliers may be replaced without much problems. However, the cluster of outliers around 17.6 s stretching from depths of 5-35 mm would prove challenging to replace. In case the identified clusters strongly correlate to turbulent structures, this filtering process might result in a bias by neglecting these structures.
In Fig. 5(a), it is also shown that there was an underestimation of the time-averaged velocity profiles by vector and sliding averaging. In order to have a closer look at this underestimation, we consider a single interrogation window at an imaging depth of approximately 12 mm (dashed line in Fig. 5(a)). In Fig. 7(d), the joint probability distribution function between the SNR of correlation planes and the "instantaneous" velocity normalized by the time-averaged velocity computed by correlation averaging is shown. For vector averaging, several outliers are present. These outliers are asymmetric about the expected value, which leads to an underestimation in time-averaged velocities. The fraction of outliers is too high to be detected and replaced by conventional techniques. For sliding averaging, the distribution of "instantaneous" velocities is a lot more favourable with a far smaller fraction of outliers. Thus, sliding averaging outperforms vector averaging.
Choosing an appropriate threshold for SNR will eliminate a fraction of the outliers. To show the effect of the choice of threshold on the time-averaged velocity computed by sliding averaging, we consider Fig. 8. Basically, we consider a relatively favourable and adverse situation of the flow. The favourable situation is at a lower velocity (0.51 m/s) and laminar conditions, while the adverse situation is at a higher velocity (1.88 m/s) and turbulent conditions.
In the more favourable scenario, the choice of the SNR threshold has virtually no effect on the time-averaged velocity profile. Moreover, there is good agreement between the velocity profiles obtained by sliding averaging and correlation averaging. As expected, the penetration depth by sliding averaging is lower than that by correlation averaging. The imposition of a SNR threshold would reject a fraction of data at each imaging depth (for example, the white dots in Fig. 7(a)-(c)). However, for the laminar flow, virtually no data are rejected at the imaging depths where sliding averaging performs well ( Fig. 8(b)).
In the more adverse scenario, the choice of the SNR threshold has a major effect on the time-averaged velocity profile. When the SNR threshold is 0, the sliding averaging results grossly underpredict the results obtained by correlation averaging, presumably because of excessive outliers. Increasing the SNR threshold to 3 improves the situation slightly. From Fig. 8(b), it can be seen that 25-60% of data are rejected. Increasing the SNR threshold to 5 greatly improves the estimation of the time-averaged velocity, and it agrees well with the results obtained by correlation averaging. However, in order to achieve this, about 80-95% of the data is rejected. As expected, the penetration depth by sliding averaging is lower than that by correlation averaging.