Turbulence statistics in a negatively buoyant multiphase plume

We investigate the turbulence statistics in a {multiphase plume made of heavy particles (particle Reynolds number at terminal velocity is 450)}. Using refractive-index-matched stereoscopic particle image velocimetry, we measure the locations of particles {whose buoyancy drives the formation of a multiphase plume,} {together with the local velocity of the induced flow in the ambient salt-water}. {Measurements in the plume centerplane exhibit self-similarity in mean flow characteristics consistent with classic integral plume theories.} The turbulence characteristics resemble those measured in a bubble plume, {including strong anisotropy in the normal Reynolds stresses. However, we observe structural differences between the two multiphase plumes. First, the skewness of the probability density function (PDF) of the axial velocity fluctuations is not that which would be predicted by simply reversing the direction of a bubble plume. Second, in contrast to a bubble plume, the particle plume has a non-negligible fluid-shear production term in the turbulent kinetic energy (TKE) budget. Third, the radial decay of all measured terms in the TKE budget is slower than those in a bubble plume.} Despite these dissimilarities, a bigger picture emerges that applies to both flows. The TKE production by particles (or bubbles) roughly balances the viscous dissipation, except near the plume centerline. The one-dimensional power-spectra of the velocity fluctuations show a -3 power-law that puts both the particle and bubble plume in a category different from single-phase shear-flow turbulence.


Introduction
Plumes containing bubbles, particles and droplets are present in both environmental and industrial applications. A few examples of environmental interest are settling sediment, volcanic eruption columns, CO 2 ocean sequestration plumes, and rising oil droplets and gas bubbles from oil well blowouts (Freeth 1987;Baines & Sparks 2005;Baines 2008;Socolofsky & Bhaumik 2008;Woods 2010;Socolofsky et al. 2011;Huppert & Neufeld 2014;Wang et al. 2016). Suspension-flow plumes differ from traditional singlefluid plumes in that the energy due to buoyant forcing is transmitted indirectly from the suspended phase to the continuous phase. The relative motion between the two phases introduces additional length and time scales; these must be included in the model formulation employed to predict plume behavior, for example, when extending classic single-phase integral plume models (Morton et al. 1956) to multiphase plumes (Milgram 1983;Sun & Faeth 1986). Identifying such scales is nontrivial due to the complexity of particle-particle and particle-fluid interactions.
For the special case of air bubbles in water, empirical data collection has allowed accurate closure of predictive schemes (Lance 1991;Risso 2018;Risso & Ellingsen 2002;Mercado et al. 2010;Almeras et al. 2017). However, the bubble-in-water plume can be quite different from other suspension-flow plumes of interest, because bubbles are far less dense than the ambient fluid (specific gravity of order 10 −3 ) and have negligible inertia. Other plumes of interest are droplets in water (specific gravity of order 1), solid particles in water (specific gravity of order 2 to 10), or liquid droplets in air (specific gravity of order 10 3 ). Each of these different suspensions can behave quite differently than the others in terms of the interphase interactions. Empirical data are relatively scarce for these other suspension-flow plumes.
Recent progress in suspension flows (especially turbulent flows) offers the hope that predictive techniques will eventually describe overall plume behavior from a direct description of the internal dynamics of particle-fluid coupling. To help support this strategy and address some of the open questions that remain about turbulent suspension flows (Balachandar & Eaton 2010;Guazzelli & Morris 2012), we take an observational approach herein, measuring the particle and fluid behaviors in a particle-laden plume. Because particles are dynamically quite different from air bubbles in water, we are curious to see how the two-phase kinematics and the interstitial fluid turbulence behave.
Both bubbles and particles have been shown to modulate the turbulent properties of multiphase flows in a manner sensitive to volume fraction (φ), particle Reynolds number (Re d ), and particle Stokes number (St). One starting point for understanding such turbulence modulation is the agitation of a quiescent fluid by bubbles or particles distributed homogeneously in space. This is often referred to as pseudo-turbulence (e.g. Lance (1991); Cartellier & Rivière (2001)) as its characteristic scales are set by the bubble or particle wakes. Pseudo-turbulence is fundamentally anisotropic, with stronger velocity fluctuations aligned in the direction of suspended-phase motion (x 1 , typically vertical). Numerical simulations of pseudo-turbulence by Riboux et al. (2013) show that the velocity fluctuation energy spectrum has a slope of k −3 for wavelengths smaller than the integral length scale; they also find that the integral length scale is the ratio of the bubble or particle diameter and the drag coefficient, C d .
A recent topical review article by Risso (2018) provides a comprehensive summary of pseudo-turbulence due to a homogeneous swarm of bubbles rising in a quiescent medium. First, the turbulence intensity is anisotropic and is dominated by the vertical component. The probability density function (PDF) of the vertical component is positively skewed due to the wake immediately behind each rising bubble, whereas the other two components are symmetrically distributed with exponential tails. Because the turbulent kinetic energy (TKE) production in suspension flows occurs at two length scales (first at the scale of bubble wakes and second at the scale of overall bubble population), the combination of the two results in a k −3 spectral subrange.
To our knowledge, there are no existing studies examining turbulence in negatively buoyant particle plumes released into an unstratified, initially quiescent fluid. However, there are many studies of turbulence in bubble plumes with notable contributions coming from Soga & Rehmann (2004); Wain & Rehmann (2005); García & García (2006); Seol et al. (2009);Simiano et al. (2009);Lai & Socolofsky (2019). These investigations reveal different velocity characteristics at different axial distance (x 1 ) from the origin. Simiano et al. (2009)  bubbles influences the effective spreading, volume fraction, and mean velocity profiles. Seol et al. (2009) and Lai & Socolofsky (2019) focused on the far field (2 < x 1 /D < 11) turbulence characteristics, such as Reynolds-stress, turbulent transport and TKE budget; their results will be used to compare with our particle plume data. Characterizing the turbulence statistics in a multiphase plume is challenging primarily due to the difficulty in simultaneously measuring velocity in both phases. Traditional intrusive techniques such as hotwire anemometry suffer from the possibility of being damaged by the solid particles. Optical measurements from techniques such as particle image velocimetry (PIV) and laser doppler anemometry are usually distorted due to the difference in optical properties of the two media. We overcome the distortion herein by carefully choosing two media with matched refractive indices. This paper is organized as follows: section 2 provides a detailed description of the plume-generating facility and method for characterizing the two-phase flow. We report and discuss our experimental observations in section 3. Finally, section 4 provides the key conclusions from this work.

Plume facility
The plume experiments were conducted in a rectangular tank (80 cm deep × 80 cm wide × 365 cm long) as described in Bordoloi & Variano (2017). The tank was filled with tap water filtered to 5 microns and maintained using an UV purifier. 91.1 kg of commercial sodium chloride (Cargill Top-Flo) was then mixed to yield a salt concentration of 0.04 g/mL. The resulting solution has density ρ s = 1.04 g/cm 3 and kinematic viscosity µ s = 1.059×10 −3 Pa.s (see Table 1).
A schematic of the experimental setup is shown in figures 1a. The negatively buoyant plume was created by releasing 110 mL of cylindrical PFA pellets (Neoflon AP-202) from a height of 56.5 cm via a screw-feed particle release mechanism (see figure 1b). The pellets are right circular cylinders with length=diameter=2 mm (see figure 1b-A). The physical properties of the solid particles and the surrounding salt solution are summarized in Table 1. Because of the hydrophobic nature of PFA, the particles tend to trap and hold air bubbles when added to water. To prevent these air bubbles from entering the experiment, particles were presoaked in water in a separate container and rapidly stirred to dislodge all bubbles. Once the particles were free of air bubbles, they were placed in a funnel for eventual release into the quiescent salt-water mixture. The funnel was kept partially submerged 23 cm below the free surface through a nozzle with internal diameter d 0 = 11.25 mm so that the particles did not contact air (see figure 1b). Particle release was governed by a motor-driven helical screw. Prior to release, the particles were held inside the funnel by the blades of the screw. Upon release, the motor rotated the screw  Figure 1: a) Schematic of experimental facility and plume release mechanism (in inset), b) regions of interest with specific dimensions, a picture of teflon particles (inset A) and an illustration of two-dimensional transformation from world coordinates into the plume coordinates (inset B).
at a constant rate of 0.5 RPM operated via Lego Mindstorms software. The particle flux, Q 0 = 7.5 cm 3 /s was measured by video-recording the release of a known quantity of particles and measuring the time difference between the exit of the first and last particle from the funnel nozzle. Before each experiment, the tank fluid was seeded with optical tracers, specifically 13 µm silver-coated hollow glass spheres (SH400S20, Conduct-O-Fil, Potters Industries). Twenty plume releases provided a sufficient number of samples for the analysis described in Section 3.
Turbulence statistics in a negatively buoyant plume

Refractive index matching
To measure the fluid velocity inside and around the plume, we matched the refractive index (n) of the particle and fluid phases. A target refractive index of PFA (n ≈ 1.34) was achieved from an aqueous solution of commercial sodium chloride with a concentration of 0.04 g/mL. The refractive index of the solution was measured to be 1.338 using an ATAGO refractometer and found to match with the salinity vs. n prediction given in Tan & Huang (2015). Although the refractive indices of the two phases were matched, because of the scattering properties of PFA the particles appeared bright under the laser illumination. To limit the intensity of particles below the saturation threshold of the camera's sensor we used a circular polarizer on each camera lens. Figure 2 shows example images of Teflon particles in different-salinity water samples illuminated by a laser sheet (wavelength=532 nm). The high-intensity sections in the images are the particles intersected by the laser sheet, and the low-intensity sections are the out-ofplane particles. Fluid tracers are much more visible when particles and fluid have matched refractive indices (figure 2b compared to figure 2a), especially behind or in front of outof-plane particles as shown in the highlighted regions.

Multiphase velocimetry
We performed stereoscopic particle image velocimetry (SPIV) to find the twodimensional three-component (2D3C) velocity field of the fluid phase of the plume. The origin of our coordinate system is located at the tank center. The plume is axisymmetric about the vertical axis, but we use Cartesian coordinates (x 1 , x 2 , x 3 ) for our measurements as described in figure 1b. Velocity measurement of the particle phase is currently underway and beyond the scope of this paper.
The x 3 = 0 plane was illuminated with a laser sheet (1 mm thick at beam waist; Quantel/Big Sky Lasers, 532-nm double-pulse Nd-YAG). Two charge-coupled device (CCD) cameras (Imager PRO-X, 1600 pixels × 1200 pixels) were synchronized with the laser pulses. They were placed ±55 • to the laser sheet (c.f. 90 • for standard 2D2C PIV). To minimize distortion due to this angle, water-filled acrylic prisms were placed between the camera lenses and the tank walls. The cameras were mounted with Nikkor Figure 3: Sample mean-subtracted instantaneous two-dimensional-three-component (2D3C) velocity field in the fluid phase superposed with the particle image in the laboratory coordinate system. The vectors show the in-plane velocity (u 1 , u 2 , m/s), and the shades of red and green colors show the out-of-plane velocity component (u 3 , m/s). The reference vector on the upper right corner corresponds to 0.5 m/s. 105 mm lenses, circular polarizers, and Scheimpflug adapters. The interframe delay (∆t) was optimized to 0.5 ms. The PIV images were collected at a frequency of 14.0 Hz. Both fluid and particle motions were decorrelated at timescales > 1/14 s. This setup gave approximately 50 independent samples during the steady-state phase of each experimental run.

Fluid-phase processing
We computed the fluid-phase velocity field using DaVis 8.2 software (LaVision GmbH). Before stereo-PIV processing, we performed stereo self-calibration to reduce disparity in the alignment between the laser sheet and the measurement plane to below 0.1 pixel. Tracers and particles in an image were separated by intensity thresholding; erosion and dilation were applied to isolate individual particles. Fluid velocity fields were obtained by masking particles and correlating tracer locations through multipass stereoscopic cross-correlation. The cross-correlation was applied with an initial interrogation window of 64 × 64 pixels and final window of 48 × 48 pixels with 50% overlap, yielding a spacing between vectors of 0.7 mm. Interrogation windows were weighted according to the symmetric two-dimensional Gaussian function. Vectors were discarded based on universal outlier detection (Westerweel & Scarano 2005) and left as data gaps. Figure 3 shows a representative field of fluid velocity fluctuations superposed with the corresponding particle image (transformed into the lab coordinate system). The grid cells are colorcoded with the out-of-plane velocity component while the vectors show the in-plane velocity field.

Particle-phase processing
To estimate the particle number density across the plume, we conducted a series of additional image-processing steps (see figure 4). First, the raw images were transformed from the camera coordinate system to the lab coordinate system using the stereoscopic Turbulence statistics in a negatively buoyant plume Figure 4: Image processing steps leading from raw PIV image to centroid identification: sample image showing a) particles from Camera 1 and Camera 2 blended into one image, b) convolution of Camera 1 image with Camera 2 image, c) image segmentation based on watershed transform, and d) final binary image with particle boundaries and located centroids in the laboratory coordinate system. mapping function. Small-size background noise features (including tracers) were removed using a simple median filter.
In stereo-mapping, after transforming into the lab coordinate system, the section of the particle intersected by the laser sheet should overlap in both cameras. Figure 4a shows a sample instant with the two camera images overlayed. The high intensity regions in the image represent the intersection between the two camera images, whereas the low intensity background is from the out-of-plane particles captured in only one of the two cameras. We remove the background noise by setting an intensity threshold and convolving the two images. The result is shown in figure 4b.
After inspecting our dataset, we find many instances where particles are nearly touching as exemplified in the lower left pair in the image. The nonuniform intensity gradients in the bordering regions among different particle pairs yielded different intensity gradients. A traditional intensity-gradient-based image segmentation therefore could not differentiate between two particles, and so adopted a segmentation technique based on the watershed transform. The watershed transform treats an image as a topographic map with brighter intensity pixels as heights, and finds the separating line that runs along the ridges (Gonzalez & Woods 2007). The result, with individual particles identified with separate colors, is shown in figure 4c. We applied an area threshold to remove partially illuminated particles and clusters that our method could not isolate. The centroids of particles identified from the final binary image are shown in figure 4d.

Coordinate transformation
Typically in jets and plumes the cross-stream velocity components (U 2 and U 3 ) are smaller than the axial velocity component (U 1 ) near the centerline. Thus, any small misalignment between the plume axis and the measurement axis could lead to a systematic bias in the measured U 2 and U 3 . The scenario is schematically shown in inset B in figure 1b in a simplified x 1 − x 2 plane in which the plume axis (x 1 ) makes an angle θ 1 with the measurement axis x 1 . For our data, if θ 1 is defined the same way and if θ 2 is the angle between the plume axis and the measurement axis on the x 1 − x 3 plane, using axisymmetry and assuming that U 2 , U 3 = 0 at the plume centerline, the measured velocities can be transformed into the plume co-ordinate system via correction angles, The measurement bias in U 2 and U 3 at the plume centerline are captured in figure 5a that shows the respective mean radial profiles, U 2 and U 3 . We perform this transformation for the 20 replicate datasets independently with correction angles (θ 1 , θ 2 ) reported in figure 5c. The bias in the non-zero centerline velocities are reflected in θ 1 and θ 2 which are small and less than 1.5 • . The two velocity components transformed into the plume coordinate systems are shown in figure 5b. The analysis that follows will be in the plume coordinate system.

Mean flow characteristics
We first characterize the mean interstitial fluid velocity based on approximately 1000 independent PIV snapshots. The mean axial velocity field U 1 (x 1 , x 2 ) in dimensional form is shown in figure 6a. The radial (x 2 ) variations of U 1 normalized by the local mean centerline velocity U c (x 1 , 0) (written as U c for simplicity from here onward) are shown in figure 6b at representative axial locations (12.9d 0 < x 1 < 15.1d 0 ) across our measurement region. The Gaussian plume radius, b g , defined as the x 2 location where U 1 = e −1 U c , is obtained via a nonlinear least-squares fit of each measured profile to a normalized Gaussian function. The data closely follow the Gaussian curve (solid line in figure 6b), which is also seen in single-and multi-phase jets and plumes (Darisse et al. 2012;Milgram 1983;Lai & Socolofsky 2019).
Next, we examine the axial decay of centerline velocity (U c ) and the axial growth of Gaussian plume radius (b g ) for a particle plume by comparing them to existing data from bubble plumes (see figure 7). For the purpose of comparison, we use two nondimensional parameters: a velocity scale (u s ) and a length scale (D = gQ 0 /4πα 2 0 u 3 s ). These are based on an integral model of bubble plumes (Bombardelli et al. 2007), in which a constant α 0 = 0.083 is used as the entrainment coefficient. The computed entrainment coefficient of our particle plume to be presented later is different from this a priori constant. The velocity scale (u s = 0.23 m/s) is the terminal velocity of a single particle in quiescent fluid and is computed using a simple drag model (Clift et al. 1978). Based on Q 0 = 7.5 cm 3 /s, the length scale D is approximately 6.8 cm. This situates our axial measurement region between 2-2.5D. The overall trends in both centerline velocity and plume radius show similarity between particle and bubble plumes. The axial variation of U c /u s in a particle plume is not measurable from our data, but we can say that it falls above the curve for bubble plumes (Lai & Socolofsky 2019) (see figure 7a). The dashed and solid lines in figure 7a are the −1/3 rd power-law fits, with A=2.05 (current data, see inset) and A=1.6 (Lai & Socolofsky (2019)). The variation of b g for all the data (including ours) is captured by the solid line (b g /D = 0.114x 1 /D) in figure 7b. Within the limited axial extent (≈ 0.5D) of our measurement, the local spreading rate (β = db g /dx 1 ) is difficult to estimate. However, justified based on the collapse of our data on the solid line in figure 7b, we will use β = 0.114 in the remaining analyses in this paper.
Profiles of mean radial ( U 2 ) and out-of-plane ( U 3 ) velocities at various axial locations are shown in figures 8a and 8b, respectively. The mean radial velocity exhibits typical jet/plume behavior: within the plume radius (|x 2 /b g | < 1) the interstitial fluid flows away from the centerline, whereas outside the plume radius (|x 2 /b g | > 1) the surrounding fluid is entrained into the plume. The mean out-of-plane velocity should be zero by design (swirling motions are not introduced by the particle feeder) but it shows some variations and asymmetry far from the plume axis. The solid line in figure 8a is a nonlinear least-squares fit to a shape function, (3.1) The shape function in equation 3.1 is obtained from the radial integration of the continuity equation written in a cylindrical coordinate system after adopting the entrainment hypothesis for jets/plumes and the Gaussian-profile assumption for U 1 (Papanicolaou & List 1988;Lee & Chu 2003). Herein, the entrainment coefficient, α=0.07 is computed as a fitting parameter in equation 3.1. The spreading vs. entrainment ratio, β/α, has been measured as 2 for pure jets (Lee & Chu 2003), and 1.2 for bubble plumes (Lai & Socolofsky 2019). Our results show β/α = 0.114/0.07 = 1.63, situating the particle-laden plume between pure jets and bubble plumes. We compute particle number density by counting the particles across each sample (detection method described in section 2.3.2). The normalized distribution of particle number density φ across the radius of the plume is shown in figure 9. Here, φ signifies the probability of finding a particle in the specified radial location, and φ c is the centerline value, which measures as 0.11 at all axial locations. The distribution of φ allows us to fit a Gaussian profile and measure the half-width (b φ ), which we use to designate the 'particlecore' from here onward. The axially averaged b φ obtained from the particle concentration profiles, when normalized by the axially averaged b g , shows that b φ /b g =0.56. This ratio for the bubble plume in Lai & Socolofsky (2019) was not measured, but in earlier studies of bubble plumes, e.g. Milgram (1983), reported values of b φ /b g are in the range of 0.8-0.9. This value is somewhat higher than our value, implying that solid particles spread less rapidly than bubbles. This may be due to the fact that rising bubbles exhibit swirling motions and experience significant lateral lift force when compared to particles (Lai & Socolofsky 2019). In addition, large-eddy simulations (LES) for a range of monodispersed and polydispersed bubble plumes in Fraga & Stoesser (2016) show that the ratio b φ /b g is very sensitive to the size distribution of bubbles across the plume. They showed that the size distribution across the plume resembles a reverse-Gaussian profile with larger bubbles populating away from the centerline. This reveals the complexity related to clearly defining a bubble-core which depends on the distribution of the number-density as well as bubble-size across the plume. A future investigation on bubble plumes examining the relationship between polydispersity and turbulence characteristics will help to better explain the above differences between the particle and bubble plumes.

Fluctuating flow characteristics
In this section, we examine the nature of the fluctuating components of velocity (u i = U i − U i ; u i,rms = u 2 i ). The normalized probability density functions (PDF) of u i at the plume centerline is shown in figure 10. The respective PDF for a bubble plume Figure 9: Normalized distribution of particle number density across the plume averaged along x 1 . The solid line is a Gaussian fit to the measured data, and the dotted line is the Gaussian fit to mean axial fluid velocity ( U 1 /U c ) from figure 6. The particle half width, b φ = 0.56b g . (Lai & Socolofsky 2019) is also included in each plot for comparison. Interestingly, the distribution of the axial velocity fluctuations is negatively skewed for a particle plume, such that axial fluctuations opposite to the direction of particle motion are more common than those along the direction of particle motion. This behavior is strikingly opposite to what is seen in bubble plumes and homogeneous bubble swarms (Riboux et al. 2010;Lai & Socolofsky 2019;Prakash et al. 2016), for which the axial velocity fluctuations moving in the direction of bubble motion are more common than those moving in the opposite direction. Another way to see this effect is that the mode of the distribution is slightly positive for the particle plume, while it is slightly negative for a bubble plume (see figure  10a). One implication of this contrasting result is that the fluid flow in a multiphase plume is sensitive to the direction of the plume with respect to gravity; particle plumes are not a simple reversal of bubble plumes.
The PDFs of radial and out-of-plane velocity fluctuations are symmetric about their means (figure 10b), and closely follow a Gaussian curve. We do not observe the prominent signatures of intermittency in the cross-stream components typically observed in bubble plumes and bubble swarms (Lai & Socolofsky 2019;Prakash et al. 2016).
In figure 11, we further examine the velocity fluctuations across the plume at different radial locations (x 2 /b g ). For this purpose, we sort the data into eight equal bins (bin width = 0.16b g ) for x 2 0 after checking that the choice of bin width does not change the results. Across the plume, the radial (u 2 ) and out-of-plane (u 3 ) velocity fluctuations are symmetric with nearly zero skewness (see figure 11b). The axial velocity fluctuations (u 1 ) switch from negatively skewed at the centerline (as discussed earlier) to positively skewed outside the half-radius (x 2 > 0.5b g ) of the plume. These characteristics are better captured in the statistical moments (rms, skewness and kurtosis) as shown in figure 12. Figure 12b shows that the change of sign in S(u 1 ) occurs at three-quarters of the plume width b g , which also coincides with the maximum u 1,rms (see figure 12a). Both u 2,rms and u 3,rms show their maximum at the centerline, and are consistently smaller than u 1,rms .   Figure 11: Standardized PDF of a) axial: u 1 /u 1,rms , b) radial: u 2 /u 2,rms , and out-ofplane: u 3 /u 3,rms at various locations across a particle plume. The colorbar indicates the various radial locations (x 2 /b g ).
The radial variation of kurtosis K(u i ) in figure 12c captures the increasing flatness of each distribution outside the plume half-radius.

Reynolds stresses
The Reynolds stresses across the plume are reported at all measured axial locations in figure 13a (normal stresses) and 13b (shear stress in the measurement plane). The turbulent kinetic energy (k) based on 13a is shown separately in figure 13c. A nonlinear least-squares fit of the data to a shape-preserving function (see Appendix A) captures well each profile (see fitted lines in figure 13). Results show that the turbulent kinetic energy is primarily dominated by the axial velocity fluctuations, and that it increases from the plume centerline to its maximum value near the edge of the plume (x 2 ≈ 0.75b g ). The maximum TKE located at x 2 /b g = 0.75 is about 44% of U 2 c and is about 1.5 times the TKE at the plume center. The in-plane shear stress, u 1 u 2 /U 2 c , follows a trend similar to that of a turbulent jet and increases from zero at the centerline to a maximum located near x 2 /b g = 0.75 (see figure 13b). The location of maximum shear is consistent with that of u 2 1 /U 2 c , suggesting that the shear stress is dominated by the axial fluctuations u 1 .
In figure 14, we compare the in-plane turbulent intensities normalized by the local mean axial velocity, ( u 2 1 / U 1 ), with two earlier studies on bubble plumes Lai & Socolofsky 2019). For this comparison, we use the shape functions (solid and dashed lines in figure 13a. The present data and the data from Duncan et al. (2009) show reasonably similar trends with unbounded growth away from the centerline as U 1 approaches zero asymptotically outside the plume. The results from Lai & Socolofsky (2019) show somewhat higher turbulent intensity inside the plume core. Also, the growth in their normalized turbulent intensities is not unbounded. Lai & Socolofsky (2019) attribute this deviation to finite U 1 outside of the plume caused by flow recirculation set up by the tank walls. We do not observe such wall effects in our experiment. The centerline turbulent intensities for the present data are 18%, 11% and 9% for u 1 , u 2 , u 3 , respectively. The axial turbulent intensity ( u 2 1 ) at the centerline for a bubble plume (Lai & Socolofsky 2019) is also approximately twice the other two normal intensities ( u 2 2,3 ), suggesting stronger anisotropy in multiphase plumes when compared to a single-phase jets/plumes in which the ratio u 2 1 / u 2 2,3 is about 1.4 (Wang & Law 2002). to the analytical expression for M (x 1 ) for pure vertical plumes (Lee & Chu 2003),

Conservation of kinematic momentum flux of the plume
Here, x 0 = −5.6d 0 is the virtual origin of the plume, and F 0 = Q 0 g is the buoyancy flux of the particles. Our results show reasonably good agreement with this model, suggesting that the particle plume obeys the scaling law M ∼ x 4/3 1 , which is consistent with buoyancy-driven plumes. The total momentum flux is primarily contributed by the momentum flux due to the mean flow M . The remaining contribution comes from M t , and it is customary to quantify the contribution by the local momentum amplification factor γ = M/ M . For our particle plume, γ is 1.2, averaged across all axial measurements. This result shows that γ for a particle plume is larger than that  for a single-phase jet/plume (γ = 1.07 − 1.09 Wang & Law (2002)) and smaller than a bubble plume (γ = 1.4 − 1.6, Lai & Socolofsky (2019)).

Velocity triple-correlation and turbulent transport
Because the PDF of axial velocity fluctuations in a particle and a bubble plume are oppositely skewed (see figure 10a), some interesting differences between the two flows can be identified in terms of the velocity triple-correlation terms. These terms contribute to the transport of turbulent kinetic energy and thus are important in the TKE budget (section 3.7). The radial profiles of the triple-correlation of in-plane velocity fluctuations    figure 16. To compare the trends we also show the respective profiles for a single-phase jet (Darisse et al. 2012) and a bubble plume (Lai & Socolofsky 2019). The profiles extracted from the two references are multiplied by a factor for visual comparison (see legends in figure 16).
Other than u 3 1 , all triple-correlation profiles for the particle plume show trends similar to a turbulent single-phase jet. The triple-correlation profiles of the first two terms ( u 3 1 /U 3 c and u 1 u 2 2 /U 3 c ) for the bubble plume in Lai & Socolofsky (2019) are significantly different from the single-phase jet and our particle plume, and they show 10-20 times larger magnitude compared to our particle plume results (see figure 16a and 16b). The bubble plume shows a positive u 3 1 near the centerline as opposed to the negative u 3 1 observed near the core of the particle plume, due to the opposite skewness of their u 1 distributions discussed earlier (see figure 10a).
Based on the triple-correlation terms computed above, we estimate the axial and radial transport of k in figure 17a and 17b, respectively. Since the TKE in the particle plume is primarily contributed by the axial stress ( u 2 1 ), these two profiles are nearly identical to those in figures 16a and 16c. These profiles are fitted to two shape functions (see appendix A) and are shown as solid lines in figures 17a and 17b, respectively. We use them for estimating the TKE budget in section 3.7.

Mean square gradients of velocity fluctuations
To test the conditions of local isotropy and local axisymmetry (George & Hussein 1991), we compute the in-plane derivatives of all three components of velocity fluctuations employing the 2 nd -order-accurate central-difference scheme. Table 2 and figure 18 show the second moments of axial (∂/∂x 1 ) and radial (∂/∂x 2 ) derivatives of all three components of velocity. These moments clearly do not satisfy the conditions of local isotropy (e.g.  The violation of local axisymmetry was also observed for bubble plumes (Lai & Socolofsky 2019), suggesting that this is a common characteristic of multiphase plume turbulence contrary to single-phase jets.
Computing the true viscous dissipation rate (stress power) requires all nine components of the velocity gradient tensor and each component has to be adequately resolved, which is not easy to achieve experimentally. Instead, we compute the pseudo-dissipation rate Figure 16: Normalized radial profiles of triple correlations (a, b) axial transport ( u 1 u 2 i ) and (c, d) radial transport ( u 2 u 2 i ) of in-plane components of turbulent kinetic energy in a particle plume at various axial locations (x 1 /d 0 ) indicated by the colorbar compared with a single-phase jet (Darisse et al. 2012) and a bubble plume (Lai & Socolofsky 2019 (3.3) Here, ν = 1.02 × 10 −6 m 2 /s is the kinematic viscosity of saline water at 23 • C. In our  To assess the accuracy in the dissipation calculation, we adopt a method based on conservation of total kinetic energy for multiphase plume in a Lagrangian framework as outlined in Lai & Socolofsky (2019). In essence, this method estimates the balance of the Figure 19: Radial profiles of non-dimensional mean dissipation rate at various axial locations (x 1 /d 0 ) indicated by the colorbar. total kinetic energy as production minus dissipation: and calculates the correction factor f c to balance the two sides. The mean kinetic energy, K ≈ U 1 2 /2 because U 1 U 2 , U 3 . This equation can be simplified in the cylindrical coordinate system (r, θ, z) by invoking steady state as follows. The left hand side of the equation is computed using the shape functions of the mean radial profiles for U 1 and k (see section 3.1 and 3.2). The first term on the right hand side is obtained by using at multiple axial locations. The last term in equation 3.4 refers to the kinetic energy production due to the particle buoyancy flux and equals CF 0 ∆x 1 , where the prefactor dη depends on λ = b φ / b g and the momentum amplification factor γ. Based on λ = 0.6 and γ = 1.2, we obtain C = 1.44. The simplified equation for a volume bounded by two axial locations x 1 and x 2 is written as: Here, I 1 , I 2 and I 3 are the axisymmetric surface integrals (2π ∞ 0 (.) rdr) of the respective terms inside the volume integral in equation 3.4 (summarized in Table 3).
We test the balance of equation 3.5 at multiple axial locations separated by 1.5b g and have found f c to vary between 2.8 and 5.1 with a mean value of 4.0. This method estimates the amount by which is underestimated by our evaluation of equation 3.6, but it is unable to correct locally along the radius of the plume. Because of this limitation in our measurement of , we use the corrected dissipation profile only for a qualitative assessment of the TKE budget discussed in the following section. Term in eq. 3.4 Full form Similarity form Value Table 3: Simplified similarity expressions and values of various integrals used in equation 3.5. The shape functions F 4 (η) and F 5 (η) are summarized in Appendix A.

Turbulent kinetic energy budget
In the cylindrical coordinate system, the time-averaged transport equation for k of an axisymmetric, steady, turbulent plume without swirl is (Kataoka & Serizawa 1989): In addition to the terms already present in single-phase jets/plumes (A: advection, P s : shear-production, T : turbulent transport, and T p : pressure transport), a multiphase plume has another source term (P b ) that represents the interfacial energy transfer at the boundaries between the dispersed phase and the surrounding fluid. Obtaining the complete TKE budget for a particle plume is challenging for multiple reasons. First, it is difficult to accurately measure experimentally the mean dissipation rate because of the under-resolved velocity gradients (as seen in section 3.6). Second, there is no existing model for pressure-velocity correlation in multiphase flow turbulence, and therefore the pressure transport term (T p ) requires an assumption. Finally, the interfacial energy transfer (P b ) is the most complex and inaccessible term, and it requires the local distribution of particles and their relative velocity with the surrounding fluid (Santarelli et al. 2016). In an ideal scenario with all other terms measured correctly, the term P b can be sought as the closing term of equation 3.6, such that P b = − (A + P s + T + T p ). We will use this method to estimate P b from our data.
We directly compute A, P s and T from our experimental data. The expressions for these in terms of similarity variables are summarized in Appendix B. Each plot contains an error-band corresponding to the 95% confidence interval in the constituent fits. Each uncertainty band is computed via the rule of uncertainty propagation described in Charonko & Prestridge (2017). The pressure transport term (T p ), not directly available from experiment, is replaced by Lumley's model T p = −2/5T , widely used for singlephase turbulence (Lumley 1978). These four terms and the mean dissipation rate ( ) multiplied by the correction factor f c = 4 are shown in figure 20a. The remaining terms of the balance equation are the interfacial energy transfer (P b ) from particles to fluid and out-of-balance experimental error (OOBE). The budget terms for a bubble plume from Lai & Socolofsky (2019) are reproduced in figure 20b for comparison. Each budget term is normalized by U 3 c /b g ; the order-of-magnitude difference in the two budgets arises because of the difference in the downstream measurement locations with respect to the source, leading to the differences in U c and b g (see figures 7a and 7b). Below we discuss the comparison of the TKE budget for the two multiphase plumes in light of previous experiments on a single-phase jet (Lai & Socolofsky 2018) and a variable-density single-phase (VDSP) jet (Charonko & Prestridge 2017). The qualitative summary of this comparison is provided in table 4.
The shape of the advection (A) profiles are very similar across all four plumes with a positive central lobe inside the plume core. The shear production term (P s ) in the particle plume increases from zero at the centerline to a dominant peak at about r/b g = 0.65 near the edge of the plume. This behavior strongly resembles a single-phase jet both in shape and relative magnitude. The P s profile in the VDSP jet also matches this trend, except that a negative P s is observed at the jet centerline. By contrast, the bubble plume shows almost negligible P s across the radius of the plume. The turbulent transport term (T ) in the particle plume exhibits a positive peak near the centerline and a negative lobe near the edge of the plume. The bubble plume and the VDSP jet are similar, but with a larger peak in the former. The single-phase jet does not match the others, i.e. near the centerline the turbulent transport is consistently observed to be negative or marginally positive (Darisse et al. 2015). : similar profile shapes : dissimilar profile shapes Table 4: Qualitative comparison summary of three TKE budget terms in our particle plume and those in a bubble plume (Lai & Socolofsky 2019), a single-phase jet (Darisse et al. 2015;Lai & Socolofsky 2018), and a variable density single-phase (VDSP) jet (Charonko & Prestridge 2017).
Because the energy production at the two-phase boundaries (P b ) is computed indirectly, the peak location in P b is sensitive to the experimental error in . Despite this limitation, it is worth noting that in both particle and bubble plumes, the sum (A + T + P s + T p ) is significantly smaller than the mean dissipation rate outside the plume core (r > 0.5b g ). This results in the TKE budget being an approximate balance between P b and .

Two-point correlation and energy spectra
We obtain the two-point statistics from the spatial autocorrelation function R 11 computed along the axial measurements of axial velocity fluctuations, u 1 . Figure 21 shows R 11 inside the plume core (x 2 /b g < 0.5) at various radial locations (shaded circles) and their meanR 11 (red circles). The horizontal axis is normalized by the particle diameter d p . We estimate the Taylor microscale (λ f ) of this flow by fitting an osculating parabola, (Pope 2000) to the first four points ofR 11 (see solid line in figure 21), yielding λ f = 1.3d p , which is of the order of the particle size. The Taylor microscale (λ f ) is related to the mean square velocity gradient in homogeneous isotropic turbulence (Pope 2000) and we use this relation to approximate Based on this method, ∂u1 ∂x1 2 (b g /U c ) 2 at x 2 = 0 is 3.2±0.6 and at x 2 = 0.5b g is 5.8 ±0.6. Within the uncertainty bounds, these values match those obtained from direct computation of gradients in section 3.6 (see figure 18a). Risso (2018) suggests that the inter-scale energy transfer in bubble-induced agitation (BIA) is quite different from the classical picture of shear-induced turbulence (SIT). As a first step toward understanding the energy cascade in particle-plume turbulence, we generate the one-dimensional energy spectrum (E 11 ) from the autocorrelation function R 11 (Pope 2000). The normalized spectra for our particle plume is compared with experimental results (Riboux et al. 2010) and DNS simulations ) in figure 22a. The result shows striking similarity in the spectra between homogeneous bubble swarms (both experiment and DNS simulation) and our particle plume, all following the earlier prediction of a κ −3 power law for wave numbers smaller than 1/λ f . Beyond κ 1 > 1/λ f , both experimental results (ours and Riboux et al. (2010)) show different slopes than the DNS results, which can be attributed to experimental uncertainty in resolving   Figure 22: Comparison of a) one-dimensional power spectra E 11 (κ 1 ) of axial velocity fluctuations and b) one-dimensional dissipation spectra in particle plume with earlier studies related to homogeneous bubble swarms (experiment: (Riboux et al. 2010), DNS simulation: ). The horizontal axis is normalized by the wave number corresponding to respective particle diameter (d p ).
velocity fluctuations below λ f . Nonetheless, this result is valuable and it suggests powerlaw universality of multiphase flow turbulence which is quite different from that for single-phase turbulent shear flows. The one-dimensional dissipation spectra D 11 (κ 1 ) = 2νκ 2 E 11 (κ 1 ) is shown in figure  22b. Quite strikingly, the peaks of dissipation for our data and bubble plume/swarms reside at κ 1 d p ≈ 0.2 − 0.4. This indicates that in both bubble-and particle-laden turbulence, production and dissipation occur near the scale of particle size. This lack of scale separation was first postulated for homogeneous bubbly flows in Lance (1991) and is observed in direct numerical simulations reported in .

Conclusions
We report an experimental study characterizing the turbulence inside a heavy particle plume descending under gravity within a salt-water solution. We measure the three components (2D3C) of the interstitial fluid velocity and the spatial distribution of particles in the central plane of the plume using refractive-index-matched stereoscopic particle image velocimetry. Below we summarize our key findings primarily in the light of the results for a bubble plume (Lai & Socolofsky 2019).
The induced liquid flow inside the plume evolves with the mean flow characteristics of a bubble plume (Lai & Socolofsky 2019), that show Gaussian mean axial velocity ( U 1 ) profile. The radial profiles of mean particle number density are also Gaussian with a halfwidth b φ equals to 0.56b g , where b g is the Gaussian half-width of the mean streamwise velocity profile in the fluid.
The turbulence inside the plume is highly anisotropic with maximum streamwise turbulence intensity measuring up to 2.7 times that of the other two components, a result consistent with that observed in a bubble plume. The PDF of the axial velocity fluctuation (u 1 ) at the centerline of the plume is skewed so that the strongest events are in the direction opposite to mean flow. This behavior is strikingly opposite to that observed previously for a bubble plume and a homogeneous swarm of bubbles. The turbulent kinetic energy (TKE) and the in-plane shear-stress peak near ≈ 0.75b g , locating the shear layer slightly outside the edge of the particle-core (b φ ). Although these quantities for a bubble plume in Lai & Socolofsky (2019) peaked at a similar location (≈ 0.7b g ), the location of the shear layer with respect to their bubble-core (b φ ) remains unknown. Another distinction between the particle and the bubble plume is observed in the fluid shear production term (P s ) in the TKE budget. While the shear production in the bubble plume is negligible across the plume compared to the other terms, we observe a distinct peak in P s at the shear layer region of a particle plume, resembling a single-phase jet.
The other TKE terms including advection (A) and the turbulent transport (T ) in the two flows exhibit similar profiles. Further, we show that the turbulent transport (T ) term follows similar qualitative profiles for both plumes and also an earlier reported variabledensity single-phase jet (Charonko & Prestridge 2017), all of which differ from a typical single-phase jet. The difference in the relative magnitude in T at the centerline between the two plumes is attributed to the difference in their third moment (skewness) of the axial velocity fluctuations.
Despite the above differences inside the core, both particle and bubble plumes show qualitatively that the TKE production by the particles (P b ) approximately balances the mean dissipation rate ( ) away from the centerline. Further, the one-dimensional spectrum in the particle plume exhibits the −3 rd power-law consistent with bubble plume and homogeneous swarms of bubbles. These two results support the notion that there is a lack of separation between the scales of production and dissipation in multiphase turbulent flows like ours.  The in-plane shear-stress is fitted into a polynomial-Gaussian function in equation A 2.
S = u z u r /U 2 c = 0.0124η + 0.0286η 3 − 0.0113η 5 exp −1.47η 2 (A 2) The axial and radial transport of the turbulent kinetic energy (k) are fitted into the polynomial-Gaussian functions in equations A 3 and A 4, respectively.

Appendix B.
Expressions for the advection (A), shear-production (P s ), and turbulent tranport (T ) terms in cylindrical co-ordinate system in the turbulent kinetic energy budget are summarized below. Each term is expressed in terms of the similarity variable η = r/b g .