# On the Modelling of Ship Wakes in S-Band SAR Images
and an Application to Ship Identification

###### Abstract

We present a novel ship wake simulation system for generating S-band Synthetic
Aperture Radar (SAR) images, and demonstrate the use of such imagery for the
classification of ships based on their wake signatures via a deep learning
approach. Ship wakes are modeled through the linear superposition of wind-
induced sea elevation and the Kelvin wakes model of a moving ship. Our SAR
imaging simulation takes into account frequency-dependent radar parameters,
i.e., the complex dielectric constant ($\varepsilon$) and the relaxation rate
($\mu$) of seawater. The former was determined through the Debye model while
the latter was estimated for S-band SAR based on pre-existing values for the
L, C, and X-bands. The results show good agreement between simulated and real
imagery upon visual inspection. The results of implementing different training
strategies are also reported, showcasing a notable improvement in accuracy of
classifier achieved by integrating real and simulated SAR images during the
training.

Index Terms—  Ship wakes, sea waves, SAR simulation, S-band, NovaSAR-1, vessel
classification, sea modelling

##  1 Introduction

Synthetic Aperture Radar (SAR) technology exhibits weather-independent
performance and enables observations in any atmospheric and illumination
conditions [1]. SAR imaging is widely used in a variety of applications,
including maritime surveillance, environmental monitoring, and disaster
response [2, 3].

Ship wake detection in SAR images is gaining popularity as the ship wake
signatures can serve as a way to estimate ship dynamics parameters and can
also be utilized for more practical applications, e.g., building a model for
ship classification based on the wakes [4, 5]. This latter application unlocks
the capability to identify ship classes from SAR images, offering the
potential for advanced maritime surveillance capabilities. However, developing
such classification models requires a huge volume of training images, which
can be challenging to obtain due to the limited availability of publicly
accessible datasets.

The majority of SAR wake simulation literature has primarily focused on L-,
C-, and X-band radar, as these are used by some of the most widely utilised
SAR satellite missions (such as TerraSAR-X, COSMO-SkyMed, ICEYE, ALOS-2,
ERS-2, and SEASAT [6, 7, 8, 9]). However, new satellite missions such as
NovaSAR-1 and NISAR-S utilise S-band instruments due to their advantages for
maritime monitoring applications. Motivated by the above, our work focuses on
simulating ship wakes specifically as imaged under a S-Band SAR system, and
then utilizes the resulted images for classification purposes in real
scenarios. Therefore, the simulation results from this work will not only
supplement the existing modelling in terms of the frequency band used but also
provide insights into the characteristics of ship wake signatures in the
S-band scenario.

Ship wake simulation under S-Band SAR system have not been previously
conducted largely due to limited information on the frequency-dependent
parameters used for the simulation, and the absence of real datasets for
verifying simulation results. Two essential frequency-dependent parameters,
i.e., the complex dielectric constant ($\varepsilon$) and relaxation rate
($\mu$) for S-Band were estimated and provided in this work. These parameters
were then used for calculating the normalized radar cross-section (NRCS) of
sea surface and hydrodynamic modulation effect respectively. The complex
dielectric constant $\varepsilon$ for S-band SAR was obtained using the Debye
approach, as reported by Meissner and Wentz in [10, 11], which is based on an
extensive experimental dataset for seawater. Meanwhile, the value of $\mu$ was
estimated from pre-existing values in the L, C, and X-bands, according to [6]
and [12].

To demonstrate the utility of the simulated SAR images, we employed the images
to train a deep-learning model, AlexNet, for classifying real ships in SAR
images acquired by NovaSAR-1. Different training strategies, including
transfer learning and mixing of simulated and limited numbers of real images
in the training dataset, have been conducted and evaluated concerning their
impact on classification performance.

In summary, our contributions are two-fold. Firstly, to the best of our
knowledge, we are the first to conduct ship wakes simulation for S-Band SAR.
Therefore, our simulation can serve as a theoretical foundation for future
work. Secondly, the training strategies outlined in our work provide a pathway
for utilizing simulated ship wake images to build a model for classifying real
ship wakes images.To demonstrate the utility of simulated SAR images, we
employed the images to train a deep-learning model, AlexNet, for classifying
real ship images acquired by NovaSAR-1.

##  2 MATERIALS AND METHODS

This work employs prior work on ship wake generation, which we have here
extended to S-band SAR; a comprehensive explanation of the previous simulator
can be found in [6] while code can be accessed through the University of
Bristol Research Data Repository at
doi.org/10.5523/bris.el0p94vgxjhi2224bx78actb4.

###  2.1 Sea Elevation Model

In this work, the sea elevation model is formulated as the superposition of
the wind-induced sea surface and the Kelvin wakes generated by a moving ship.
Wind-induced sea elevation at time t can be expressed as [6, 12, 13]:

$\displaystyle Z_{sea}$ $\displaystyle=\sum_{i}\sum_{j}A_{ij}$
$\displaystyle\cos\left({k_{i}\left[x\cos{\theta_{i}}+y\sin{\theta_{j}}+\omega
t+\sigma_{ij}\right]}\right),$ (1)

where k and $\omega$ respectively are the wavenumber and radial frequency
while $\sigma$ is the initial phase uniformly distributed between (0, 2$\pi$).
$A_{ij}$ represents the amplitude expressed as:

$A_{ij}=\sqrt{2S(k_{i})D(k_{i}\theta_{j})\Delta{k_{i}\Delta{\theta_{j}}}},$
(2)

where $S(k)$ and $D(k_{i}\theta_{j})$ are wave spectrum and angular spreading
function [6].

In modelling the waves generated by moving ships, we used the Kelvin waves
model,$Z_{ship}$, represented by Zilman et al. in [14]:

$Z_{ship}=\frac{V_{s}}{g}\frac{\partial{\Phi_{s}}}{\partial{x}},$ (3)

with $V_{s}$ and $\Phi_{s}$ being ship velocity and an approximated form of
fluid velocity potential. Finally, composite sea elevation model ($Z$) due to
existing wind and moving ship can be obtained as the superposition of
$Z_{sea}$ and $Z_{ship}$.

###  2.2 Frequency-Dependent Parameters

There are two essential parameters that need to be determined to model the
interaction between incoming waves and the sea surface, i.e., complex
dielectric constant ($\varepsilon$) and the relaxation rate ($\mu$).

The complex dielectric constant ($\varepsilon$) has been calculated based on
the work of Meissner and Wentz in [10, 11]. Therein, extensive experiments
were conducted to determine $\varepsilon$ for both pure and sea water over a
temperature range from -20°C to +40°C and a frequency range up to 90 GHz. The
dielectric constant ($\varepsilon$) of sea water with a specific salinity
($S$) and temperature ($T$), radiated under a specific wavelength ($\lambda$),
can be calculated as:

$\displaystyle\varepsilon(\lambda,S,T)$ $\displaystyle=$
$\displaystyle\varepsilon_{R}+\frac{\varepsilon_{S}(S,T)-\varepsilon_{R}}{1+{\left[i\lambda_{R}(S,T)/\lambda\right]}^{1-\alpha}}-2\frac{i\sigma(S,T)\lambda}{c},$
(4)

Here, $\lambda_{R}$, $\varepsilon_{S}$, and $\sigma$ are all dependent on $S$
and $T$, while $\varepsilon_{R}$ and $\alpha$ are set to 4.44 and 0.012,
respectively, and $c$ represents the speed of light. The detailed
implementation of above equation can be found on [10]. Assuming that $S=35$
ppt [15], $T=21^{\circ}C$ [16], and NovaSAR-1 frequency is 3.2 GHz [17], the
corresponding dielectric constant for our simulation is 69.63-i38.95.

The magnitude of $\mu$ is poorly known experimentally, especially for S-Band
case. In our experiment, the $\mu$ value for S-band was interpolated from the
known $\mu$ values for L, C, and X-band, which were provided in [12, 18].
Finally, the estimated value for $\mu$ in the NovaSAR-1 band are $\mu=0.05$
s-1 for $V_{w}$ $\leq$ 5 m/s and $\mu=0.39$ s -1 for $V_{w}$ $>$ 5 m/s.

###  2.3 SAR Imaging of Sea Surface

To simulate SAR images of the modeled sea surface, we employ the two-scale
model (TSM). This model is based on the resonant Bragg scattering theory
describing the scattering of electromagnetic waves from a rough surface [6,
13]. The Bragg scattering solution for normalized radar cross-section (NRCS)
with VV and HH polarizations is expressed as [6]:

$\sigma_{0}=8\pi k_{e}^{4}\cos^{4}{\theta_{l}}W(k_{Bx},k_{By})|T_{c}^{2}|,$
(5)

where $k_{e}$ is radar wavenumber. The $\theta_{l}$, $W$, and $T_{c}$ are
respectively local incidence angle, wavenumber spectral density and complex
scattering function. Then, the NRCS after taking into account hydrodynamic and
tilt modulations:

$\sigma(x,y)=\sigma_{0}\left[1+\int{(M(k)F(k)e^{ikx}+c.c)dk}\right],$ (6)

where $F(k)$ is the 2D FFT of composite sea elevation model $Z$. $M(k)$ is
defined as follow [12, 19]:

$M(k)=-4.5\frac{k_{y}^{2}}{|k|}\frac{\omega-i\mu}{\omega+i\mu}+\frac{4\cot{\theta_{r}}}{1\pm\sin^{2}{\theta_{r}}}ik_{y},$
(7)

where $\theta_{r}$ is the nominal radar incident angle. The SAR image
intensity can then be obtained using:

$\displaystyle I(x,y)$ $\displaystyle=\iint{\partial(y_{i}-y)}$
$\displaystyle\frac{\sigma(x,y)}{P_{a}^{{}^{\prime}}(x,y)}{\exp{\left(-\pi^{2}{\left[{\frac{x_{i}-x-\frac{R}{V}U_{r}}{b}}\right]}^{2}\right)}}dxdy.$
(8)

$\partial(.)$ is the Dirac delta function of the impulse response function in
range direction while $R$, $H$, and $V$ are the range distance, platform
height, and platform velocity, respectively. $P_{a}^{{}^{\prime}}$ is degraded
azimuth resolution which depends on number of incoherent looks and radar
wavelength [6]. Finally, multiplicative random noise is added to the final
image to simulate speckle, using the Weibull probability density function
(PDF) as described in [20].

###  2.4 Ship Classification

In this work, the AlexNet [21] model was trained to differentiate wakes
produced by cargo and tanker ships. Due to the limited availability of real
ship wake images for testing, we utilized ship wake crops from five real
NovaSAR-1 products and applied data augmentation by rotating, flipping, and
varying noise amplitude. In total, this dataset consisted of 960 real images
and is hence denoted as NVReal dataset.

Three different strategies for training the model have been conducted. Both
Experiment I and Experiment II are based on a transfer learning approach. For
this purpose, a pre-trained model was produced using a large dataset,
SynthwakeSAR, a synthetic dataset contains a total of 46,080 simulated SAR
ship wakes images [22]. SynthwakeSAR dataset is accessible at
doi.org/10.5523/bris.30kvuvmatwzij2mz1573zqumfx. Since the SynthwakeSAR
dataset was originally generated for X-band SAR, we reproduced the dataset to
represent images acquired by an S-band SAR system (denoted as SynthwakeSAR-S
dataset in this work) and used that for initial training. Afterwards, a
simulated NovaSAR-1 dataset (NVSim) was generated based on parameters of real
ship wake images from NovaSAR-1 (NVReal). This NVSim dataset was used to re-
train the pre-trained model via fine-tuning (FT) for Experiment I and via
feature extraction (FE) schemas for Experiment II.

NVSim contains 18,000 images which were produced by assigning parameters
within the range of the real parameters for the real NovaSAR images in NVReal,
as shown in the Table 1 (denoted as R superscripts). This step was to ensure
that the simulated images closely resemble their real counterparts so that the
classification model can learn their representation for recognizing the real
images. Individual ship parameters are provided in Table 2.

Table 1: Sea state and imaging parameters for generating NVSim. Parameters |  Symbol |  Value range |  Number of sampling
---|---|---|---
Wind speed 10m | $U_{10}$ |  $U_{10}^{R}\pm 0.5$ m/s | 3
Wind direction | $\theta_{w}$ | $\theta_{w}^{R}\pm 15^{\circ}$ | 10
Ship speed | $V_{s}$ |  $V_{s}^{R}\pm 0.5$ m/s | 5
Ship direction | $\theta_{s}$ | $\theta_{s}^{R}\pm 30^{\circ}$ | 8
Incident angle | $\theta_{i}$ | $\theta_{i}^{R}\pm 2^{\circ}$ | 3
Table 2: Ship parameters for generating NVSim. MMSI |  Class |  $L(m)$ |  $B(m)$ |  $T(m)$ |  $Fetch(km)$ |  $U_{10}^{R}$ (m/s) |  $\theta_{w}^{R}$(∘) |  $V_{s}^{R}$ (m/s) |  $\theta_{s}^{R}$(∘) |  $\theta_{i}^{R}$(∘)
---|---|---|---|---|---|---|---|---|---|---
215071000 | Tanker | 183 | 32 | 9.30 | 80 | 4.13 | 343.50 | 6.73 | 351.50 | 29.35
218433000 | Cargo | 126 | 23 | 6.90 | 25 | 6.10 | 359.50 | 7.45 | 182.50 | 29.31
371720000 | Tanker | 229 | 103 | 9.50 | 80 | 3.33 | 332.50 | 7.56 | 355.50 | 23.96
636013817 | Tanker | 227 | 150 | 9.70 | 80 | 4.85 | 4.85 | 6.89 | 156.50 | 27.87
477665200 | Cargo | 266 | 150 | 10.50 | 80 | 6.20 | 341.50 | 8.44 | 181.50 | 28.23

In contrast to Experiment I and Experiment II, Experiment III did not employ
any pre-trained network. In this case, the model was directly trained using
the NVSim dataset. This experiment aims to assess the model’s performance in
the absence of a pre-existing dataset. Furthermore, we incorporated a portion
of images from NVReal and mixed them with NVSim during the training stage for
all experiments. The amount of real images involved was varied from 10% to
30%. We did not extend beyond this range as it would reduce the number of the
remaining NVReal for testing the model. The flowchart of the experiments is
given in Figure 1.

Fig. 1: Flowchart of experiments. The SynthwakeSAR-S refers to reproduced
version of SynthwakeSAR using NovaSAR-1 imaging properties while NVSim and
NVReal are respectively simulated and real datasets of ship wakes.

##  3 RESULTS AND DISCUSSION

In this section, we present the visual comparison between the real and
simulated images. Moreover, the performance of the classification model
trained under different strategies is also discussed.

###  3.1 Simulated NovaSAR-1 Images

In general, evaluating the accuracy of simulated images can be achieved
through both qualitative and quantitative approaches. However, for our case,
like-for-like comparison between simulated and real images is challenging
because simulated images may not be entirely identical but rather akin or
resembling the real ones. Therefore, measuring similarity between the images,
for example using Structural Similarity (SSIM) index or correlation
coefficient, would be inappropriate as these metrics work on a pixel-by-pixel
basis.

Considering the aforementioned limitation, we chose to compare the images
through visual inspection. As depicted in Figure 2, the simulated image
closely resembles the real one, except for the appearance of turbulent wake on
the real images. Furthermore, as the ship moves faster and parallel to the
platform direction, the height of wind-induced ambient waves increases,
resulting in the transverse waves becoming scarcely visible. This phenomenon
leaves only the prominent Kelvin arms visible in both images, indicating a
consistent pattern or agreement between the two images. This behavior has also
been noted in [6, 23].

Fig. 2: Comparison between simulated (A) and real (B) images of ship wakes for
NovaSAR-1 SAR system. The parameters used are: $L$=126 m, $B$=23 m, $T$=6.9 m,
$\theta_{s}=182^{\circ}$, $V_{s}=7.45$ m/s, $\theta_{wind}=359.5^{\circ}$,
$V_{w}=6.10$ m/s, and $\theta_{i}=29.31^{\circ}$.

###  3.2 Classification Results and Performance

In this section, the use of simulated images for classifying real images is
demonstrated. Table 3 shows the classification accuracy for each experimental
setup under different proportions of NVReal included for training the model.

Table 3: Classification performance for all three experiments. NVReal Portion (%) | Accuracy (%)
---|---
| Exp. I | Exp. II | Exp. III
0 | 60.41 | 50.86 | 51.56
10 | 60.25 | 60.41 | 63.96
20 | 40.02 | 64.14 | 59.97
30 | 68.55 | 72.13 | 71.83

Based on the experimental results, three key findings have been highlighted.
Firstly, the model accuracy is generally low if no real images are introduced
to the model during training. While it may be assumed that this result
suggests a poor representation of real images in the simulated dataset, it is
noteworthy that the accuracy across all experiments remains above 50%. This
implies that training the model solely using NVSim does provide actionable
information and is not reduced to random guessing.

Secondly, introducing a portion of real images in the training process can
significantly boost the model performance. It is found that adding 30% of the
available NVReal data to NVSim in the training process improves the accuracy
by almost 12% for Experiment II. This result indicates that it is still
possible to build a good classification model even when there is limited real
data available for training.

Lastly, fine-tuning a model pre-trained under synthetic ship wake dataset
(i.e, Experiment II) is capable of producing a higher accuracy of up to 72%.
This result suggests that combining representative simulated dataset and re-
training all the network weights could provide benefits for building a
classification model in the absence of a real dataset for training.

It should be noted that our results are based on experiments conducted on five
ship models due to the challenges in obtaining real images along with their
ship and sea states, as well as imaging parameters. Therefore, a more refined
result is likely to be obtained with an increased number of real datasets.
This would enhance the model’s ability to learn features from a wider range of
real images rather than relying extensively on augmented images. Nevertheless,
our results are presented to showcase the benefits of using simulated images
in developing a deep learning model for classifying real ships in SAR images.

##  4 CONCLUSIONS

In this work, the first results on simulation of ship wake acquired under
S-Band SAR system is reported. The simulation has been conducted by employing
a two-scale model, which requires frequency-dependent parameters, i.e., the
complex dielectric constant ($\varepsilon$) and the relaxation rate ($\mu$) of
seawater. The value of $\varepsilon$ was determined using the Debye model,
whereas $\mu$ was estimated for the S-band based on pre-existing values in the
L, C, and X-bands.

Based on visual comparison, there is good agreement between the simulated and
real images. Additionally, we showcased the utilization of simulated images to
address a real-world challenge, i.e., ship wake classification on NovaSAR-1
images. We have found that a classification model having accuracy of up to 72%
can be attained by incorporating a very small proportion of real images along
with simulated ones in the training stage.

## References

  * [1] Y. S. Tsai, A. Dietz, N. Oppelt, and C. Kuenzer,  “Remote Sensing of Snow Cover Using Spaceborne SAR: A Review,”  Remote Sensing, vol. 11, no. 12, 2019.
  * [2] S. Adeli, B. Salehi, M. Mahdianpari, L. J. Quackenbush, B. Brisco, H. Tamiminia, and S. Shaw,  “Wetland Monitoring Using SAR Data: A Meta-Analysis and Comprehensive Review,”  Remote Sensing, vol. 12, no. 14, 2020.
  * [3] E. Makhoul, S. V. Baumgartner, M. Jäger, and A. Broquetas,  “Multichannel SAR-GMTI in Maritime Scenarios With F-SAR and TerraSAR-X Sensors,”  IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens., vol. 8, no. 11, pp. 5052–5067, 2015.
  * [4] B. Tings, S. Wiehle, and S. Jacobsen,  “Ship Wake Component Detectability on Synthetic Aperture Radar (SAR),”  in IGARSS 2020 - 2020 IEEE Int. Geosci. and Remote Sens. Symp., 2020, pp. 1233–1235.
  * [5] P. Heiselberg, K. Sørensen, and H. Heiselberg,  “Ship Velocity Estimation in SAR Images Using Multitask Deep Learning,”  Remote Sens. Environ., vol. 288, pp. 113492, 2023.
  * [6] I. G. Rizaev, O. Karakuş, S. J. Hogan, and A. Achim,  “Modeling and SAR Imaging of the Sea Surface: A Review of the State-of-the-art with Simulations,”  ISPRS J. Photogramm. Remote Sens., vol. 187, pp. 120–140, 2022.
  * [7] J.K.E. Tunaley, E.H. Buller, K.H. Wu, and M.T. Rey,  “The Simulation of the SAR Image of a Ship Wake,”  IEEE Trans. Geosci. Remote Sens., vol. 29, no. 1, pp. 149–156, 1991.
  * [8] F. Nunziata, A. Gambardella, and M. Migliaccio,  “An Educational SAR Sea Surface Waves Simulator,”  Int. J. Remote Sens., vol. 29, no. 11, pp. 3051–3066, 2008.
  * [9] I. Hennings, R. Romeiser, W. Alpers and A. Viola,  “Radar Imaging of Kelvin Arms of Ship Wakes,”  Int. J. Remote Sens., vol. 20, no. 13, pp. 2519–2543, 1999.
  * [10] T. Meissner and F.J. Wentz,  “Algorithm Theoretical Basis Document (ATBD) Version 2 :AMSR Ocean Algorithm ,”  EOS Project Goddard Space Flight Center, NASA, pp. 1–3, 2000.
  * [11] T. Meissner and F.J. Wentz,  “The Complex Dielectric Constant of Pure and Sea Water from Microwave Satellite Observations,”  IEEE Trans. Geosci. Remote Sens., vol. 42, no. 9, pp. 1836–1849, 2004.
  * [12] C. Brüning, R. Schmidt, and W. Alpers,  “Estimation of the Ocean Wave-radar Modulation Transfer Function from Synthetic Aperture Radar Imagery,”  J. Geophys. Res., vol. 99, no. C5, pp. 9803–9815, May 1994.
  * [13] W. Ren, P. Liu, X. Ren, and Y. Jin,  “SAR Image Simulation of Ship Turbulent Wake Using Semi-Empirical Energy Spectrum,”  IEEE J. Multiscale Multiphys. Comput. Tech., vol. 6, pp. 1–7, 2021.
  * [14] G. Zilman, A. Zapolski, and M. Marom,  “On Detectability of a Ship’s Kelvin Wake in Simulated SAR Images of Rough Sea Surface,”  IEEE Trans. Geosci. Remote Sens., vol. 53, no. 2, pp. 609–619, 2015.
  * [15] W. J. Gould and S. A. Cunningham,  “Global-scale Patterns of Observed Sea Surface Salinity Intensified Since the 1870s,”  Commun. Earth Environ., vol. 2, no. 1, pp. 76, Apr 2021.
  * [16] “Global Sea Surface Temperature Reaches a Record High,” https://climate.copernicus.eu/global-sea-surface-temperature-reaches-record-high, 2023,  [Online; accessed 25-December-2023].
  * [17] R. Bird, P. Whittaker, B. Stern, N. Angli, M. Cohen, and R. Guida,  “NovaSAR-S: A Low Cost Approach to SAR Applications,”  in Proc. of 2013 Asia-Pacific Conf. on Synthetic Aperture Radar (APSAR), 2013, pp. 84–87.
  * [18] E. A. Caponi, D. R. Crawford, H. C. Yuen, and P. G. Saffman,  “Modulation of Radar Backscatter from the Ocean by a Variable Surface Current,”  J. Geophys. Res., vol. 93, no. C10, pp. 12,249–12,263, Oct 1988.
  * [19] D. R. Lyzenga,  “Numerical Simulation of Synthetic Aperture Radar Image Spectra for Ocean Waves,”  IEEE Trans. Geosci. Remote Sens., vol. GE-24, no. 6, pp. 863–872, 1986.
  * [20] O. Karakuş, E. E. Kuruoğlu, and A. Achim,  “A Generalized Gaussian Extension to the Rician Distribution for SAR Image Modeling,”  IEEE Trans. Geosci. Remote Sens., vol. 60, pp. 1–15, 2022.
  * [21] A. Krizhevsky, I. Sutskever, and G. E. Hinton,  “ImageNet Classification with Deep Convolutional Neural Networks,”  in Adv. Neural Inf. Process. Syst. 2012, vol. 25, pp. 1097–1105, Curran Associates, Inc.
  * [22] I. G. Rizaev and A. Achim,  “SynthWakeSAR: A Synthetic SAR Dataset for Deep Learning Classification of Ships at Sea,”  Remote Sensing, vol. 14, no. 16, 2022.
  * [23] B. Tings,  “Non-Linear Modeling of Detectability of Ship Wake Components in Dependency to Influencing Parameters Using Spaceborne X-Band SAR,”  Remote Sensing, vol. 13, no. 2, 2021.


# Effects of surfactants on bubble-induced turbulence

Tian Ma1<EMAIL_ADDRESS>Hendrik Hessenkemper1<EMAIL_ADDRESS>Dirk
Lucas1  Andrew D. Bragg2<EMAIL_ADDRESS>1Helmholtz-Zentrum Dresden –
Rossendorf, Institute of Fluid Dynamics, 01328 Dresden, Germany 2Department of
Civil and Environmental Engineering, Duke University, Durham, NC 27708, USA

###### Abstract

We use experiments to explore the effect of surfactants on bubble-induced
turbulence (BIT) at different scales, considering how the bubbles affect the
flow kinetic energy, anisotropy and extreme events. To this end, high-
resolution Particle Shadow Velocimetry measurements are carried out in a
bubble column in which the flow is generated by bubble swarms rising in water
for two different bubble diameters ($3$ mm $\&$ $4$ mm) and moderate gas
volume fractions ($0.5\%\sim 1.3\%$). To contaminate the flow, different
amounts of 1-Pentanol were added to the flow, leading to different bubble
shapes and surface boundary conditions. The results reveal that with
increasing surfactant concentration, the BIT generated increases in strength,
even though bubbles of a given size rise more slowly with surfactants. We also
find that the level of anisotropy in the flow is enhanced with increasing
surfactant concentration for bubbles of the same size, and that for the same
surfactant concentration, smaller bubbles generate stronger anisotropy in the
flow. Concerning the intermittency quantified by the normalized probability
density functions of the fluid velocity increments, our results indicate that
extreme values in the velocity increments become more probable with decreasing
surfactant concentration for cases with smaller bubbles and low gas void
fraction, while the effect of the surfactant is much weaker for cases with
larger bubble and higher void fractions.

###### keywords:

##  1 Introduction

Surfactants are ‘surface-active’ molecules and/or particles that are easily
adsorbed at surfaces and form monolayers (Langevin, 2014; Manikantan &
Squires, 2020). They are present in most multiphase systems, being either
naturally present or else introduced purposefully (Levich, 1962; Stone, 1994;
Soligo et al., 2019; Lohse, 2022). In bubbly flows, it has been shown that
even small amounts of surfactant can cause dramatic changes to the bubble
shapes (Tomiyama et al., 2002; Tagawa et al., 2014; Hessenkemper et al.,
2021a); bubble rise velocities (Bel Fdhila & Duineveld, 1996; Cuenot et al.,
1997; Takagi et al., 2008; Tagawa et al., 2014), lateral migration (Lu et al.,
2017; Hayashi & Tomiyama, 2018; Ahmed et al., 2020; Atasi et al., 2022),
cluster formation (Takagi et al., 2009; Maeda et al., 2021), coalescence
(Verschoof et al., 2016; Lohse, 2018; Néel & Deike, 2021) and mass transfer on
the interfaces (Cuenot et al., 1997; Schlüter et al., 2021).

The mechanism by which surfactants influence the velocity field in the
vicinity of a gas-liquid interface was first introduced by Frumkin & Levich
(1947) and Levich (1962), who showed that as a bubble rises in an aqueous
surfactant solution, surfactant is swept off the front part of the bubble by
surface convection and accumulates in the rear region. They lower the surface
tension in the rear region relative to that at the front, leading to a
gradient of surface tension along the interface. This gradient creates a
tangential shear stress on the bubble surface (Marangoni effect) which opposes
the surface flow, causes the interface to become more rigid, and increases the
drag coefficient $C_{D}$. The free-slip boundary condition that occurs at a
gas-liquid interface for an ideal purified system (e.g. ‘hyper clean’ water)
breaks down, and the rise speed of the contaminated bubble decreases with
increasing surfactant concentration, approaching the behaviour of a rigid body
for sufficient surfactant concentration.

The effect of drag coefficient enhancement due to the adsorption of
surfactants has been studied in great detail. Readers are referred to
Magnaudet & Eames (2000), Palaparthi et al. (2006) and Takagi & Matsumoto
(2011) for detailed reviews from the perspective of hydrodynamics, and
Manikantan & Squires (2020) from the perspective of surfactant dynamics. An
early study (Bachhuber & Sanford, 1974) measured the terminal velocity of
small bubbles (bubble Reynolds number $Re_{p}\sim O(10)$) rising in distilled
water at two heights in the flow and observed that the velocity at the lower
height was in good agreement with that expected for a clean spherical bubble,
while that measured at the greater height was consistent with a particle
having a value of $C_{D}$ that corresponds to a rigid sphere. Their
interesting observation reflects at least two well-known aspects of how
surfactants impact bubble motion. First, water used in typical lab conditions
is contaminated, possibly containing considerable amounts of surfactants that
can influence the motion of bubbles. Second, it takes a finite time (that
depends on the properties of the surfactants) for the surfactants to be
adsorbed at the bubble surface before reaching an equilibrium state, with the
implication that bubble rise velocities decrease with increasing travel
distance until reaching a constant value (Durst et al., 1986; Tagawa et al.,
2014; Hessenkemper et al., 2021a).

The fact that the bubble rise velocity decreases in the presence of
surfactants also leads to a smaller inertial force experienced by bubbles. As
a result of this, for a fixed bubble size, the bubble shape is less flattened
in a contaminated system than it is in a purified system, despite the fact
that surface tension is reduced by surfactants.

The aforementioned modifications in bubble rise velocity, surface boundary
condition and shape due to surfactants have a strong impact on the wake
structure and the path instability of a rising bubble. Tagawa et al. (2014)
used experiments to investigate the effect of surfactants on the path
instability of an air bubble rising in quiescent water, and categorized the
rising bubble trajectories as straight, helical, or zigzag based on the bubble
Reynolds number and surface slip condition. They observed that the trajectory
of the bubble was first helical and then transitioned to a zigzag motion in
Triton X-100 solution – a phenomenon never reported in purified water (Mougin
& Magnaudet, 2001; Shew & Pinton, 2006). Pesci et al. (2018) used numerical
simulations to study the effects of soluble surfactants on the dynamics of a
single bubble rising in a large spherical domain. They found that if the
surface contamination is sufficiently high, a quasi-steady bubble velocity can
be obtained over a wide range of surfactant concentrations, independent of the
exact concentration value in the bulk. Furthermore, they also observed a
transition from a helical to a zig-zag rising bubble trajectory as reported
experimentally by Tagawa et al. (2014), but also found that the nature of the
trajectory depends on both the initial surface and bulk surfactant
contamination. Pesci et al. also looked at the vorticity in the vicinity of
the bubble and observed strong vorticity production very close to the bubble
surface due to Marangoni forces, while in clean water this behaviour does not
appear for path unstable bubbles (Mougin & Magnaudet, 2006). Legendre et al.
(2009) performed numerical simulations to study the two-dimensional flow past
a cylinder and investigated the influence of a generic slip boundary condition
on the wake dynamics. They showed that slip markedly decreases the vorticity
intensity in the wake. Mclaughlin (1996) and Cuenot et al. (1997) used a
stagnant-cap approximation to compare the wake structure produced by
contaminated bubbles and solid spheres. The former study revealed that the
wake volume for contaminated bubbles is larger than for solid spheres moving
at the same Reynolds number, and the latter study found that the wake length
is larger for contaminated bubbles. They explained that the increase of the
wake size is caused by the abrupt change to the dynamic boundary condition
where the transition from a quasi-shear-free (the upper part of the bubble) to
a quasi-no-slip (the rear region) condition generates strong vorticity at the
interface. As a result, there is more vorticity injected into the flow for
contaminated bubbles than for a rigid sphere with a uniform no-slip condition,
resulting in a larger wake.

The observations summarized above were mainly for isolated bubbles in the
absence of turbulence. Swarms of bubbles can however induce strong background
turbulence, and the observations above could indicate that turbulence
generated by wakes in bubble swarms could depend strongly on the degree of
contamination in the fluid. Direct numerical simulations (DNS) have revealed
the following properties of dilute dispersed bubbly flows: (i) the liquid
velocity fluctuations are highly anisotropic, with fluctuations that are much
larger in the direction of the mean bubble motion (Lu & Tryggvason, 2013; Ma
et al., 2020; du Cluzeau et al., 2022). It was recently demonstrated that this
strong anisotropy also exists at the small-scales of BIT due to the energy
being injected at the scale of the bubble (Ma et al., 2021). (ii) The
probability density functions (PDFs) of all fluctuating velocity components
are non-Gaussian, with the PDF of the vertical velocity fluctuations being
strongly positively skewed, while the other two directions have symmetric PDFs
(Roghair et al., 2011; Riboux et al., 2013). (iii) There is strong enhancement
of the turbulent kinetic energy dissipation rate in the vicinity of the bubble
surface (Santarelli et al., 2016; du Cluzeau et al., 2019). (iv) The energy
spectra in either the frequency or wave number space exhibit a $-3$ power law
scaling for a subrange for all the components of the fluctuating fluid
velocity (Pandey et al., 2020; Innocenti et al., 2021). (v) Pandey et al.
(2020) and Innocenti et al. (2021) showed that on average the energy transfer
is from large to small scales in BIT, however, there is also evidence of an
upscale transfer when considering the transfer of energy associated with
particular components of the velocity field (Ma et al., 2021). (vi) For a
bubbly flow with background turbulence, the flow intermittency is
significantly increased by the addition of bubbles to the flow when compared
to the corresponding unladen turbulent flow with similar bulk Reynolds number
(Biferale et al., 2012; Ma et al., 2021). Note that some of these DNS studies
effectively considered contaminated bubbly flows since they used no-slip
conditions on rigid bubble surfaces, while others used slip boundary
conditions, and a clear understanding of exactly how the aforementioned
properties of BIT depend on contaminants in the flow is not available. For
further details on the relevant DNS studies, readers are referred to the
recent reviews of Risso (2018) and Mathai et al. (2020).

Several experiments have also observed these properties for bubbly turbulent
flows (Lance & Bataille, 1991; Rensen et al., 2005; Riboux et al., 2010;
Mendez-Diaz et al., 2013; Lai & Socolofsky, 2019; Masuk et al., 2021).
However, in most of these experiments there will be some contamination in the
liquid, and the impact of this on the results is not understood. Indeed, a
systematic investigation into how contamination affects the properties of BIT
is still missing, to the best of the authors’ knowledge. Therefore, in the
present paper we systematically explore the effect of surfactants on the
properties of BIT produced by bubble swarms, considering both single-point and
two-point turbulence statistics to characterize large and small scale flow
properties. The rest of this paper is organized as follows. In § 2, we
describe the experimental set-up and the measurement techniques. We then give
a brief overview regarding the effect of the surfactant on the single bubble
behaviour in the chosen solution (§ 3), followed by a presentation of the
single-point statistics for the bubble swarm in § 4. Finally, the multipoint
results from the experiments that give insights into the properties of the
flow at different scales are divided into two parts, namely, the flow
anisotropy in § 5.1 and the extreme events in § 5.2.

##  2 Experimental method

###  2.1 Experimental facility

Figure 1: Sketch of the bubble column used in the experiments (note that in
the actual experiment, the number of bubbles in the column is $O(10^{3})$).
The sketch is not to scale; the column depth is many times larger than the
bubble diameter. Right shows an instantaneous realization of velocity vector
over the FOV in the case LaTap with two in-focus bubbles recognizable by their
sharp interfaces and associated wakes.

The experimental apparatus is identical to that in Ma et al. (2022), and we
therefore refer the reader to that paper for additional details; here, we
summarize. Figure 1 shows a sketch of the experimental set-up, consisting of a
rectangular column (depth $50\,\mathrm{mm}$ and width $112.5\,\mathrm{mm}$)
made of acrylic glass with a water fill height of $1,100\,\mathrm{mm}$. Air
bubbles are injected through $11$ spargers which are homogeneously distributed
at the bottom of the column.

We use tap water in the present work as the base liquid and add 1-Pentanol as
an additional surfactant with varying bulk concentration $C_{\infty}$ of $0$
ppm, $333$ ppm and $1,000$ ppm. The surfactants properties and the single-
bubble behaviour in these solutions will be discussed in detail in § 3.2. Note
that the tap water will already be sightly contaminated prior to adding the
surfactants, and the bubbles can behave differently in this tap water without
surfactants compared to that in a pure water system (no surfactants, see e.g.
Veldhuis et al., 2008; Takagi & Matsumoto, 2011). For the bubble sizes
considered in the present study (with $d_{p}>2$ mm), however, the slight
contamination in the tap water does not have a significant effect on the
bubble motion (Ellingsen & Risso, 2001). This point was confirmed as well by
our previous study (Hessenkemper et al., 2021b) for the bubble sizes
considered here, with both the bubble rise velocity and bubble aspect ratio
showing little difference between tap water and purified water systems.

We consider two different bubble sizes by using spargers with different inner
diameters. For each bubble size, we maintain the same gas inlet velocity and
ensure that all cases are not in the heterogeneous regime of dispersed bubbly
flows. In total, we have six mono-dispersed cases labelled as SmTap, SmPen,
SmPen+, LaTap, LaPen, and LaPen+ in table 1, including some basic
characteristic dimensionless numbers for the bubbles. Here, “Sm/La” stand for
smaller/larger bubbles and “Tap/Pen/Pen+” stand for corresponding cases,
having 1-Pentanol concentration with $C_{\infty}=$ $0$ ppm, $333$ ppm and
$1,000$ ppm, respectively. It should be noted that the three cases with larger
bubble sizes have higher gas void fraction than the three cases with smaller
bubbles. This is because in our setup it is not possible to have the same flow
rate for two different spargers while also maintaining a homogeneous gas
distribution for mono-dispersed bubbles.

Parameter | SmTap | SmPen | SmPen+ | LaTap | LaPen | LaPen+
---|---|---|---|---|---|---
$C_{\infty}$ (ppm) | 0 | 333 | 1000 | 0 | 333 | 1000
$\alpha$ | 0.46% | 0.47% | 0.54% | 1.36% | 1.33% | 1.33%
$d_{p}\;(\mathrm{mm})$ | 3 | 2.5 | 2.6 | 4.1 | 3.8 | 3.8
$\chi$ | 1.8 | 1.3 | 1.2 | 1.8 | 1.5 | 1.3
$Ga$ | 504 | 394 | 420 | 815 | 729 | 719
$Eo$ | 1.27 | 0.91 | 0.99 | 2.4 | 2.07 | 2.03
$Re_{p}$ | 798 | 587 | 528 | 986 | 866 | 817
$C_{D}$ | 0.531 | 0.599 | 0.841 | 0.900 | 0.933 | 1.018
Table 1: Selected characteristics of the six bubble swarm cases investigated.
Here, $C_{\infty}$ is the bulk concentration of 1-Pentanol, $\alpha$ the
averaged gas void fraction, $d_{p}$ the equivalent bubble diameter, $\chi$ the
aspect ratio, $Ga\equiv\sqrt{\left|\pi_{\rho}-1\right|gd_{p}^{3}}/\nu$ the
Galileo number, $Eo\equiv\Delta\rho gd_{p}^{2}/\sigma$, the Eötvös number. The
values of $Re_{p}$, the bubble Reynolds number and $C_{D}$, the drag
coefficient are based on $d_{p}$ and the bubble to fluid relative velocity
obtained from the experiment.

###  2.2 Flow imaging

#### Liquid phase

For all measurements stated below, we use a $2.5$ megapixel CMOS camera
(Imaging Solutions) equipped with a $100$ mm focal length macro lense
(Samyang). To measure the liquid velocity we use Particle Shadow Velocimetry
(PSV), which is similar to planar Particle Image Velocimetry (PIV) with the
only difference being that backlight illumination together with a shallow
Depth of Field (DoF) is used. The velocity is then determined by correlating
the displacement of sharp tracer particles inside a narrow sharpness region. A
detailed description of the image processing procedure used can be found in
Hessenkemper & Ziegenhein (2018). The measurement setup and data acquisition
are similar to our previous work (Ma et al., 2022), so that only the key
aspects for the liquid velocity measurements are stated in what follows.

The liquid velocity measurements take place along the $x_{1}-x_{2}$ symmetry
plane in the centre of the depth ($x_{3}$). The measurement height with
$x_{1}=0.65$ m is based on the centre of the field of view (FOV) (figure 1).
We use $10\,\mathrm{\mu m}$ hollow glass spheres (Dantec) as tracer particles
with estimated Stokes number $O(10^{-3})$, and they are hydrophilic and so are
not adsorbed by the bubble surface. The flow section is illuminated with a
$200$ W LED-lamp and an f-stop of $2.8$ is set at the camera lens to provide
an effective shallow DoF $\sim 370\,\mathrm{\mu m}$. The captured images cover
a FOV of $18.8\,\mathrm{mm}\,(H_{2})\times 12.5\,\mathrm{mm}\,(H_{1})$, which
results in a pixel size of $9.8\,\mathrm{\mu m}$. For each case, $15000$ image
pairs are recorded, with a time delay of about $0.5$ sec before the next image
pair is acquired, providing us $15000$ uncorrelated velocity fields. The
comparatively large bubble shadows in the images are masked and the
corresponding areas are not considered in the following correlation step. The
final interrogation window is $64\times 64$ pixels with $50\%$ overlap,
resulting in a vector spacing of $0.627$ mm. Here, we use standard PIV-
processing steps to determine the velocity, including multi-pass/window
refinement steps and universal outlier detection. A representative transient
FOV for the case LaTap is shown in figure 1, overlaid with the resulting
liquid velocity vector field. In this figure, two in-focus bubbles can be
identified by their sharp interfaces and the associated wakes in the velocity
vector field.

#### Gas phase

Figure 2: (a) Example image of the bubbles from a fragment in the middle of
figure 3(b), identified by the area marked with white dashed line there for
the LaTap case. (b \- d) Three steps to reconstruct hidden bubble parts for an
irregular shaped bubble: (b) segmentation mask, (c) radial distances, and (d)
ground truth radial distances.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Example images of the bubbles with fitted ellipses for an arbitrary
instant. Left: smaller bubbles with (a) Case SmTap, (c) Case SmPen and (e)
Case SmPen+. Right: larger bubbles with (b) Case LaTap, (d) Case LaPen and (f)
Case LaPen+.

Similar to our previous work (Ma et al., 2022), compared to the liquid phase
we use a larger FOV $(40\,\mathrm{mm}\times 80\,\mathrm{mm})$ for the gas
phase to capture the bubble statistics in all the cases considered. For
sufficient statistics, $1500$ image pairs are evaluated at a frame rate of
$500$ f.p.s.. The bubbles are detected with a novel convolutional neural
network (CNN) based bubble identification algorithm, which is trained to
segment bubbles in crowded situations (coloured patches in figure 2b).
Furthermore, the hidden part of a partly occluded bubble is estimated with an
additionally trained neural network. For this, an equal number of radial rays
are generated for each identified segment, starting from the segment centre to
its boundary (green and red lines in figure 2c). Afterwards, radial rays that
touch neighbouring segments (red lines in figure 2c) are corrected with the
neural network in order to account for the occluded part of the bubble. These
steps are illustrated in figure 2, and we refer the reader to the paper
(Hessenkemper et al., 2022) for more technical detail. Here, we show in figure
3 the example images of detected bubbles with corrected outline for all the
six cases. While the $SmTap$, $LaTap$ and $LaPen$ cases have mostly irregular
bubble shapes associating with shape oscillations (see supplementary movie),
the cases $SmPen$, $SmPen+$ and $LaPen+$ have fixed shapes. The bubble size is
calculated using the volume-equivalent bubble diameter of a spheroid as
$d_{p}=(d_{maj}^{2}d_{min})^{1/3}$, where $d_{maj}$ and $d_{min}$ are the
lengths of the major and minor axes of the fitted ellipse, respectively. We
observe that the surfactant changes the bubble shape dramatically. For the
smaller bubbles with $C_{\infty}=1000$, the shape is close to a sphere with a
small aspect ratio $\chi=d_{maj}/d_{min}=1.2$, while for tap water case, i.e.
$C_{\infty}=0$ we obtain $\chi=1.8$. For the larger bubbles we also have
aspect ratios in a similar range, with $\chi$ decreasing from $1.8$ to $1.3$
when going from higher to lower surfactant concentrations.

After the bubble detection, the centroids are tracked in each image pair to
obtain corresponding bubble velocities. Again a shallow DoF is used, which
allows to select only sharp bubbles in the column centre by evaluating the
grey value derivative along their contour. We choose a grey value derivative
threshold that provides a centre region of about $15$ mm thickness
$(3d_{p}\sim 4d_{p})$ in which the bubbles are evaluated. Some of the
important bubble properties are listed in table 1.

It should be noted here that in the case of ellipsoidal bubbles, i.e. the tap
water cases, errors in the determined size are larger due to the bubbles being
tilted with respect to the camera. This error can result in a bubble size
overprediction of up to $25\%$ for a strongly tilted spheroidal bubble with
fixed shape, which is, however, on average much lower due to the normal
distribution of the orientation angle with a maximum around zero (Bröder &
Sommerfeld, 2007). Using the single bubble data with two cameras (separate
experiments reported in § 3.2), we estimate the overprediction of the average
volume-equivalent diameter for a fixed spheroidal shaped bubble (corresponding
to Case SmTap) to be about $4\%$. Due to the irregular (wobbling) bubble
surfaces, the error for the larger bubbles in the LaTap case are unknown, but
we expect it to be of the same order of magnitude. The error for the void
fraction is about $5-8\%$ (Hessenkemper et al., 2022) and the error of the
bubble velocity is around $4\%$.

##  3 Surfactant properties and single bubble behaviour

###  3.1 Surfactant properties

The effect of surfactants on the bubbles arises fundamentally due to their
impact on the shear stress at the gas-liquid interface, and depends upon the
species of surfactants along with their concentration. The variety of
surfactants is huge and our quantitative understanding of their effect on the
base-fluid is still in its infancy. There are also sometimes “untypical”
scenarios reported. Ybert & di Meglio (1998) showed the effect of surfactant
desorption, leading to a remobilization of the interface by considering
bubbles contaminated with sodium dodecyl sulfate (SDS). This surfactant has a
significant desorption velocity, and they found that due to the progressive
remobilization of the interface, the rise velocity of the SDS-preloaded bubble
increases over a large portion of its trajectory (see their figure 12), until
reaching a constant value. Similar phenomena was also reported by the same
group (Ybert & di Meglio, 2000) for short-chain alcohols, where they
illustrated bubble rise velocities in ethanol-water solution that were almost
indistinguishable from those in pure water, because of its fast desorption
kinetics.

In the present study, we focus on the limit of high bubble Péclet number, $Pe$
($=d_{p}\left\|\boldsymbol{U}^{G}\right\|/D$, where $d_{p}$ is the bubble
diameter, $\boldsymbol{U}^{G}$ is the averaged bubble velocity and $D$ is the
diffusion coefficient of the surfactant in the liquid), relatively low
surfactant concentration compared to the critical micelle concentration (CMC)
and low rate of desorption. Under these conditions, surface convection is
dominant compared to the adsorption–desorption kinetics and diffusive
transport of surfactants on the bubble surface. Surfactants collect in a
stagnant cap at the back end of the bubble while the front end is stress free
and mobile. Palaparthi et al. (2006) call this the “stagnant cap regime” –
that most commonly realized in typical cases of bubbles moving in a
contaminated solution, and, in this regime the bubbles are significantly
affected by the surfactant, as shown by numerical and experimental results of
the rise velocity and wake structure.

In our experiments we use 1-Pentanol as the surfactant and use bulk
concentrations $C_{\infty}$ of 0 ppm, 333 ppm and 1000 ppm for the six bubble
swarm cases. Tests for a single bubble rising in the column were also
conducted with concentrations up to 2000 ppm (discussed in § 3.2). Here, the
main reason for choosing 1-Pentanol is that with this surfactant (and these
concentrations) the bubbles adapt very quickly to the surfactants and no
transitions in their motion type appear within the FOV (Tagawa et al., 2014).

We measure the surface tension separately with a bubble pressure tensiometer,
allowing a dynamic determination of the surface tension. For 1-Pentanol
concentration of 1000 ppm, using profile analysis tensiometry (Eftekhari et
al., 2021) we observe a surface tension reduction $\sim 7\%$ compared to tap
water. The Péclet number is of order of $10^{5}$, indicting the convection
time scale is very small compared to the time scale of the diffusion on the
bubble surface. The solubility for 1-Pentanol is 2500 $\mathrm{mol/m}^{3}$,
which is much higher than the present cases (1000 ppm $\sim
9.19\,\mathrm{mol/m}^{3}$). Micelles do not form in the bulk aqueous phase as
indicted by our experimental results to be reported later for a single bubble.

###  3.2 Preliminary test for single bubble

To provide reference cases later for the bubble swarms analysis, we first
examine the effect of the surfactants on an isolated bubble in the same setup
by varying the concentration of 1-Pentanol $C_{\infty}$ with levels 0 ppm, 333
ppm, 666 ppm, 1000 ppm, 1500 ppm and 2000 ppm in tap water. This is a wider
range than that used for the bubble swarms to check if any interface
remobilization occurs at high concentration. We use a single sparger in the
column centre, but vary two types of spargers sizes for each $C_{\infty}$,
which are also used for the bubble swarm results to generate (approximately)
two bubble sizes. Here, we continuously generate single bubbles by setting a
small gas flow rate as we also did in our previous investigation (Hessenkemper
et al., 2021a). The low generation frequency of 1 Hz ensures a large enough
distance between successive bubbles to avoid any influence from a leading
bubble, but allows to study multiple same sized bubbles under the same
conditions for better statistics.

For the evaluation of the single bubble rising trajectory in three-dimensions
(3D), we conducted stereoscopic measurements with an additional second camera
of the same type placed perpendicular with respect to the imaging direction of
the first camera. The single bubble images are also captured with a frame rate
of $500$ f.p.s., but with a larger f-stop of $11$, since sharp bubble outlines
at all depth positions are required. As no overlapping bubbles have to be
distinguished in the single bubble experiments, we use a conventional image
processing approach here instead of the CNN-based one that is used for bubble
swarms. At first, sharp edges marking the outline of the bubbles are detected
with a Canny-edge detector. The solid of revolution is then used to determine
the volume of the bubble. Further geometric properties like the bubble major
axis, defined as the largest extent in a $2$D projection, and the bubble minor
axis, defined as the largest extent perpendicular to the bubble major axis,
are extracted as well from the bubble projections (Ziegenhein & Lucas, 2017).
In comparison to the usual approach of fitting ellipses around the contour of
detected bubbles, almost the same volume and semi-axis length are obtained,
with only minor deviations for irregular shaped bubbles (Ziegenhein & Lucas,
2019). Afterwards, the centroid of the bubble projections is tracked through
successive images to obtain time-resolved instantaneous bubble velocities
$\tilde{\boldsymbol{u}}^{G}$, which can be decomposed into a ensemble-averaged
part $\boldsymbol{U}^{G}$ and a fluctuating part $\boldsymbol{u}^{G}$. The
same decomposition is also used for the liquid phase with
$\tilde{\boldsymbol{u}}^{L}=\boldsymbol{U}^{L}+\boldsymbol{u}^{L}$ in the
later sections.

(a)

(b)

(c)

Figure 4: Measured single bubble size (a), aspect ratio (b) and rising
velocity (c) in different 1-Pentanol concentration.

Figures 4(a) shows the bubble diameters using two types of spargers size as a
function of $C_{\infty}$. The bubble size is slightly reduced when adding
1-Pentanol for both types of sparger, generating smaller bubbles from 2.9 mm
($C_{\infty}=0$ ppm) to 2.7 mm ($C_{\infty}=1000$ ppm) and larger bubbles from
4.4 mm ($C_{\infty}=0$ ppm) to 4.1 mm ($C_{\infty}=1000$ ppm). This is due to
the influence of the surfactants that reduce the surface tension and hence
affect the bubble formation at the rigid orifice (Loubière & Hébrard, 2004;
Drenckhan & Saint-Jalmes, 2015). A similar trend can also be found in table 1
for the corresponding bubble swarm cases. In contrast to bubble size, the
bubble shape changes dramatically (figure 4b), with the aspect ratio reducing
from $1.8$ at $0$ ppm to $1.15$ at $1000$ ppm for the smaller bubble and
$1.9\,(0\,\mathrm{ppm})\rightarrow 1.36\,(1000\,\mathrm{ppm})$ for the larger
bubble. Further increasing $C_{\infty}$ does not change the bubble shapes for
the present smaller and larger bubbles. A similar trend can be found for the
averaged rise velocity $U^{G}_{1}$ for the both bubble sizes, namely,
increasing 1-Pentanol concentration above 1000 ppm does not affect
$U^{G}_{1}$. For both bubble sizes and for all concentrations considered we
found no increase of the rise velocity, nor any change in terms of aspect
ratio above $C_{\infty}=1000$ ppm, hence, no surface remobilization has
occurred. Based on these we assume that for both bubble sizes a saturated
contamination state (at least from the hydrodynamic perspective) is reached at
the threshold $C_{\infty}\approx 1000$ ppm.

We now look more in detail at the single bubble cases with $C_{\infty}=0$ ppm,
$C_{\infty}=333$ ppm and $C_{\infty}=1000$ ppm, since they have the same set
of $C_{\infty}$ values as the swarm cases we will discuss later. We label
these cases S-SmTap, S-SmPen and S-SmPen+ for the smaller bubble and S-LaTap,
S-LaPen and S-LaPen+ for the larger bubble, respectively (the first letter S
denotes a single bubble). Figure 5(a) shows a 3D view of the bubble
trajectories for the three smaller bubble cases, while figure 5(b) shows a
view of these trajectories from the top. The trajectories are similar to those
found by Tagawa et al. (2014) for 2 mm bubbles, with helical trajectories at
lower 1-Pentanol concentration (see S-SmTap and S-SmPen) and zigzag rising
paths for higher concentration (S-SmPen+). These zigzag and helical motions
are accompanied by oscillations in the vertical velocity $\tilde{u}^{G}_{1}$
plot of figure 5(c) (corresponding to the paths in figure 5a), i.e. the
bubbles alternately speed up and slow down as they rise. As expected, the
S-SmPen+ case has much higher oscillation in $\tilde{u}^{G}_{1}$, compared to
the S-SmTap and S-SmPen cases. This is caused by the relatively unstable wake
structure that is associated with zigzag paths, rather than the more stable
one associated with the helical trajectories (Cano-Lozano et al., 2016).
Moreover, we observe the helical trajectories of S-SmTap and S-SmPen are not
as regular as those in Tagawa et al. (2014) due to our larger bubble size
$\sim 3$ mm. This is also indicated by the vertical velocities (figure 5c),
showing oscillations in S-SmTap and S-SmPen compared to the constant
$\tilde{u}^{G}_{1}$ for the cases with relatively perfect helical rising paths
in Tagawa et al. (2014).

(a)

(b)

(c)

(d)

(e)

Figure 5: 3D trajectories (a), top view of those (b) and corresponding
instantaneous bubble velocities with components $\tilde{u}^{G}_{1}$ (c),
$\tilde{u}^{G}_{2}$ (d), $\tilde{u}^{G}_{3}$ (e) over time (normalized by
$\Delta t=2$ ms) of smaller single bubbles at $C_{\infty}=0$ ppm (S-SmTap),
$C_{\infty}=333$ ppm (S-SmPen) and $C_{\infty}=1000$ ppm (S-SmPen+). All
figures are from the same track of the particular case.

(a)

(b)

(c)

(d)

(e)

Figure 6: 3D trajectories (a), top view of those (b) and corresponding
instantaneous bubble velocities with components $\tilde{u}^{G}_{1}$ (c),
$\tilde{u}^{G}_{2}$ (d), $\tilde{u}^{G}_{3}$ (e) over time (normalized by
$\Delta t=2$ ms) of larger single bubbles at $C_{\infty}=0$ ppm (S-LaTap),
$C_{\infty}=333$ ppm (S-LaPen) and $C_{\infty}=1000$ ppm (S-LaPen+). All
figures are from the same track of the particular case.

Another important observation from figures 5(c-e) is that the frequency of the
vertical velocity is approximately twice as large as that of the horizontal
velocities due to the frequency of the force oscillations in the corresponding
directions (Mougin & Magnaudet, 2006). Moreover, there is a trend that with
increasing $C_{\infty}$, the frequency of the oscillation in velocities
increases. This is in line with the observation reported in Tagawa et al.
(2014) for their smaller bubbles.

In figure 6, we plot the same quantities for the three larger bubble cases. A
fixed path type cannot be found for Case S-LaTap. Although in figure 6(e) it
shows a zigzagging trend, there are many other snapshots, showing a flattened
helix (not shown here). Case S-LaTap belongs to the chaotic regime due to the
very large $Re_{p}$ and $Eo$, while S-LaPen exhibits a flattened helical
motion, and S-LaPen+ converges toward a zigzag path. Moreover, compared to the
smaller bubbles, the vertical velocities (figure 6a) in all three larger
bubble cases display irregular oscillations.

##  4 Flow characterization for bubble swarm

(a)

(b)

Figure 7: Gas void fraction (a) and liquid/gas vertical velocity (b) along the
horizontal axis of FOV.

Basic statistics of both phases for the considered cases are plotted in figure
7. Both the void fraction (figure 7a) and the vertical gas/liquid velocity
(figure 7b) have very flat profiles, indicating statistical homogeneity of the
current flow in the FOV. Even with zero-bulk-flow (averaged over the entire
flow cross section), we observe for almost all the cases $U^{L}_{1}\neq 0$
along the horizontal axis of the FOV (especially for the cases SmPen+ and
LaPen+), so that the relative velocity in the FOV is in general not equal to
the bubble terminal rise velocity. Furthermore, due to the Marangoni effect we
find that the bubble swarms rise more slowly with increasing $C_{\infty}$ for
both smaller and larger bubbles, consistent with our results for the
corresponding single bubble in § 3.2. It is worth noting although LaPen and
LaPen+ have a similar $U^{G}_{1}$, LaPen+ has a smaller relative velocity due
to the higher $U^{L}_{1}$, as indicated by its smaller bubble Reynolds number
(table 1).

We now turn to consider the role played by the surfactant in generating BIT.
Hereafter, all average quantities refer to the liquid, so that the upper index
$L$ is dropped for simplicity. In figure 8(a) we plot the turbulent kinetic
energy (TKE), $k$, calculated as

$k=\frac{1}{2}\left((u_{1}^{rms})^{2}+2(u_{2}^{rms})^{2}\right)\;,$ (1)

assuming the out-of-plane velocity variance is equal to the measured
horizontal component. This approximation – axisymmetry about the vertical
direction is expected for the BIT dominated flows far from the wall. Here,
$u_{1}^{rms}$ and $u_{2}^{rms}$ are the root-mean-square values of the
vertical and the horizontal velocity fluctuation, respectively. For both
smaller and larger bubbles, surprisingly, the TKE is highest for the cases
with the highest $C_{\infty}$, although SmPen+ and LaPen+ have the lowest
$Re_{p}$ in the cases of smaller and larger bubbles, respectively. It is also
interesting to note that TKE in the cases SmPen and LaPen do not change much
compared to SmTap and LaTap, respectively, indicating that the amount of
1-Pentanol added ($C_{\infty}=333$ ppm) is not enough to modify the BIT
already initiated in the tap water system. However, we should keep in mind
that the bubble sizes in SmPen and LaPen are slightly smaller than their
corresponding tap water cases (see table 1), so that if the TKE is almost the
same for the $C_{\infty}=0$ ppm and $C_{\infty}=333$ ppm cases, then the
surfactant is in fact leading to a positive contribution to the BIT
generation.

In table 2 we summarise the effect of the surfactant on three aspects that
influence BIT. Aside from reducing the bubble Reynolds number, the bubble
surface instability (e.g. deformation and wobbling) also reduces when
increasing $C_{\infty}$, and both of these reductions have a negative impact
on BIT production. This suggests that the change of boundary condition induced
by increasing $C_{\infty}$ plays the key role in causing the TKE to be
enhanced by the addition of surfactants.

Factor influencing BIT | Change by increasing $C_{\infty}$ | Contribution to BIT
---|---|---
Bubble Reynolds number | reduced | negative
Boundary condition | free-slip $\rightarrow$ no-slip | positive
Surface instability | reduced | negative
Table 2: Surfactant effect on multiple factors that influence bubble-induced
turbulence.

Following our previous study (Ma et al., 2021), we define a Reynolds number
$Re_{H_{2}}\equiv u^{\ast}H_{2}/\nu$, indicating the range of scales in the
turbulent bubbly flows. Here, $u^{\ast}\equiv\sqrt{(2/3)k_{\mathrm{FOV}}}$,
and $k_{\mathrm{FOV}}$ is the TKE averaged over the FOV of the liquid phase.
In figure 8(b) we depict $Re_{H_{2}}$ versus the large-scale anisotropy ratio
$u_{1}^{rms}/u_{2}^{rms}$ (also averaged over the FOV). The figure reflects
the same behaviour as the TKE, namely the case with the highest $C_{\infty}$
also has the largest $Re_{H_{2}}$ for the corresponding bubble size group.
Furthermore, we find that for large scales whose fluctuating velocities are
characterized by $u_{1}^{rms}$ and $u_{2}^{rms}$, the smaller bubbles produce
more anisotropy in the flow than the larger bubbles, as reflected by a larger
ratio of $u_{1}^{rms}/u_{2}^{rms}$ for the cases SmTap, SmPen and SmPen+. This
is in very close agreement with our previous study based on DNS data of
bubble-laden turbulent channel flow driven by a vertical pressure gradient and
with no-slip boundary conditions on the bubble surfaces (Ma et al., 2021).

The more important point prompted by figure 8(b) is that the large-scale
anisotropy increases significantly with increasing surfactant concentration
for both smaller and larger bubbles. This behaviour is due to the effect of
the surfactant on both the wake structure and trajectory type of the bubbles
as $C_{\infty}$ changes, as was already seen for different single bubbles in
figures 5 and 6.

(a)

(b)

Figure 8: (a) Turbulent kinetic energy along the horizontal axis of FOV. (b)
Reynolds number, $Re_{H_{2}}$ plotted versus large-scale anisotropy ratio,
$u_{1}^{rms}/u_{2}^{rms}$.

##  5 Turbulence modification across scales

###  5.1 Turbulence anisotropy

The components of the second-order velocity structure function are defined as

$D_{2}^{ij}(\boldsymbol{x},\boldsymbol{r},t)\equiv\langle\Delta
u_{i}(\boldsymbol{x},\boldsymbol{r},t)\Delta
u_{j}(\boldsymbol{x},\boldsymbol{r},t)\rangle\;,$ (2)

where $\Delta u_{i}(\boldsymbol{x},\boldsymbol{r},t)$ denotes the difference
in the velocity at positions $\boldsymbol{x}$ and
$\boldsymbol{x}+\boldsymbol{r}$ at time $t$. Hereafter, we suppress the space
$\boldsymbol{x}$ and time $t$ arguments since we are considering a flow which
is statistically homogeneous and stationary over the FOV. The PSV measurement
provides access to data associated with separations along two directions,
namely, the vertical separation $\boldsymbol{r}=r_{1}\boldsymbol{e}_{1}$
($r_{1}=\|\boldsymbol{r}\|$) and the horizontal separation
$\boldsymbol{r}=r_{2}\boldsymbol{e}_{2}$ ($r_{2}=\|\boldsymbol{r}\|$). Hence,
we are able to compute the four contributions
$D^{L}_{2}(r_{1})=D_{2}^{11}(r_{1})$, $D^{L}_{2}(r_{2})=D_{2}^{22}(r_{2})$,
$D^{T}_{2}(r_{1})=D_{2}^{22}(r_{1})$, and $D^{T}_{2}(r_{2})=D_{2}^{11}(r_{2})$
based on the Cartesian coordinate system depicted in figure 1.

Figures 9(a) and (b) show the measured transverse and longitudinal second-
order structure functions of the $u_{1}$ component, respectively. The results
show that the values of the structure functions increase in the order SmTap,
SmPen, SmPen+, LaTap, LaPen, LaPen+, which corresponds to larger bubble size
and higher surfactant concentration. This ordering also holds for the $u_{2}$
component computed (not shown). Similar to the results for the TKE, while the
difference between Sm(La)Tap and Sm(La)Pen is small, $D_{2}^{\gamma\gamma}$
(no index summation is implied for $\gamma$) for Sm(La)Pen+ have noticeably
higher values across most scales for both smaller and larger bubble sizes.

(a)

(b)

Figure 9: Transverse (a) and longitudinal (b) second-order structure function
of the $u_{1}$ component, with separations along the horizontal (a) and the
vertical (b) directions.

To characterizing the multiscale anisotropy associated with $D_{2}^{ij}$, we
plot the ratios of components $D^{L}_{2}(r_{1})/D^{L}_{2}(r_{2})$ and
$D^{T}_{2}(r_{2})/D^{T}_{2}(r_{1})$ in figure 10, which would be equal to
unity for an isotropic flow (Carter & Coletti, 2017). Generally, the results
show that both ratios decrease monotonically for decreasing separation.
Moreover, the ratio of the transverse structure functions departs more
strongly from unity than the longitudinal ones. By comparing the cases with
the same bubble size, we find that generally the cases with higher surfactant
concentration generate stronger anisotropy in the flow across the scales, and
this trend is most obvious in the plot for
$D^{T}_{2}(r_{2})/D^{T}_{2}(r_{1})$. This is in agreement with the results of
the large-scale anisotropy in § 4.

(a)

(b)

Figure 10: Ratio of longitudinal (a) and transverse (b) structure functions in
different separation directions for all the cases. In (a,b) the horizontal
lines indicate the value of unity.

Furthermore, considering the cases with the same $C_{\infty}$, our results
indicate that smaller bubbles generate stronger anisotropy in the flow,
consistent with our previous studies (Ma et al., 2021, 2022) that only
considered fully-contaminated bubbles. However, the results in figure 10 show
that it is not true in general that smaller bubbles generate stronger
anisotropy in the flow, because depends on the surfactant concentration. For
example, figure 10 shows that $LaPen+$ can be more anisotropic than $SmTap$.

###  5.2 Extreme fluctuations in the flow

(a)

(b)

Figure 11: Snapshot of the original velocity vector (a) and the intensity
distributions of the normalized velocity increment $\left|\Delta
u_{1}(r_{2}=\Delta)\right|/\sigma_{\Delta u_{1}}$ (b) from the same instant
based on the SmPen+ case. The in-focus bubbles are denoted in (a).

Having explored the role of surfactants on the flow anisotropy, we now turn to
consider the effect of surfactants on extreme fluctuations of the velocity
increments – a phenomenon associated with internal intermittency in single-
phase turbulence (Frisch, 1995; Sreenivasan & Antonia, 1997). In the present
bubbly flows, extreme events in the liquid phase result from the boundary of
the wakes produced by the bubbles (Ma et al., 2022). These contributions can
be seen in figure 11, with the normalized velocity increments showing large
values at the edge of the wakes, while the velocities are largest inside the
wakes. Another possible contribution to extreme events would come from the
flow in the boundary layer at the top of the rising bubble where the flow may
abruptly change from being equal to the bubble surface velocity to the
background bulk velocity over a thin boundary layer. For spherical bubbles
with large $Re_{p}$ rising in a quiescent flow the thickness of this boundary
layer is $O(d_{p}Re_{p}^{-1/2})$ (Moore, 1963; Batchelor, 1967). However, we
cannot resolve the velocity fluctuations in this thin boundary layer with our
current experimental method (and the results discussed below should be
interpreted accordingly), and further work is needed to understand how these
regions might contribute to extreme fluctuations in the flow.

In figure 12 we plot the PDFs of the velocity increments for separations equal
to ten PSV grids $(r=10\Delta=0.18H_{2})$, corresponding to an ‘eddy size’
with $O(d_{p})$. (The results for other separations show qualitatively similar
trends and so are not shown). First, all of the six cases show that the
velocity increments have strongly non-Gaussian PDFs, just as in single-phase
turbulence. Second, while the PDFs of the transverse velocity increments in
figure 12(c,d) are almost symmetric, the PDFs of the longitudinal velocity
increments in plots (a,b) are negatively skewed for the horizontal separations
$r_{2}$. Finally, with respect to the effect of bubble size and surfactant
concentration, the results show very interesting behaviours: for the smaller
bubbles, both the longitudinal and the transverse PDFs (figure 12a,c) become
increasingly non-Gaussian in the order of SmPen+, SmPen, SmTap, which
corresponds to the order of decreasing $C_{\infty}$ of 1-Pentanol. This is
consistent with our previous finding (Ma et al., 2022) that the non-
Gaussianity of the PDFs of the velocity increments becomes stronger as the
Reynolds number (here, $Re_{H_{2}}$) is decreased in bubble-laden turbulent
flow considered here, while the opposite occurs for single-phase turbulence
(Frisch, 1995). An explanation for why the intermittency is largest for the
cases with lower surfactants concentration can be given as follows. For these
three cases the volume fraction is not large, and there are relatively few
regions in the flow where turbulence is produced due to the bubble wakes,
meaning that turbulence in the flow is very patchy and therefore intermittent.
For the cases with higher $C_{\infty}$, they have larger flow Reynolds numbers
accompanied by longer and wider bubble wakes. As a result the wake regions are
more space filling and hence the flow is less intermittent than the cases with
lower surfactant concentration which have shorter and narrower bubble wakes.
Again, this difference compared to intermittency in single-phase turbulence
can be traced back to the fact that the regions of intense small-scale
velocity increments occur near the turbulent/non-turbulent interface at the
boundary of the bubble wake, which is more similar to so-called external
intermittency (Townsend, 1949) than internal intermittency.

(a)

(b)

(c)

(d)

Figure 12: Normalized probability density functions of the longitudinal (a,b)
and transverse (c,d) velocity increments for the three smaller bubble cases
(a,c) and three larger bubble cases (b,d) with the separations along the
horizontal direction ($r_{2}=10\Delta$, where $\Delta$ is one PSV grid).

The $C_{\infty}$-dependency of the small-scale intermittency just discussed
for the smaller bubbles is however much weaker for the larger bubbles (figure
12b,d). A potential reason for this difference is that the three larger bubble
cases have volume fractions that are almost three times larger than those for
the smaller bubbles, and the bubble wake volume is also much larger, with much
higher turbulence intensity generated by the wakes. As a result of these
properties, while there are significsant regions of the flow that are almost
quiescent for the smaller bubble cases the same is not true for the larger
bubbles, with turbulent activity filling a significant fraction of the flow.
Due to this, the velocity difference in and outside the wake is smaller for
the larger bubble cases than the smaller bubble cases, and the intermittency
is less dependent on $C_{\infty}$ since even for the largest $C_{\infty}$
case, these velocity differences across the wak boundaries are already
smaller. This then explains why LaTap, LaPen and LaPen+ do not show the same
$C_{\infty}$-dependence in the PDFs as the smaller bubble cases.

While the PDFs results present the effect of the surfactants on the extreme
events at a fixed separation, we look at now this property across all the
scales quantified by the flatness $D^{T}_{4}(r_{2})/(D^{T}_{2}(r_{2}))^{2}$
and $D^{L}_{4}(r_{2})/(D^{L}_{2}(r_{2}))^{2}$ (figure 13). The first
observation from the results in figure 13 is that all six bubble-laden cases
show a similar behaviour as $r$ increases, with values of up to 30 gradually
reducing as $r$ increases. However, there is still considerable deviation from
the Gaussian value of 3 even at the largest scale $r_{2}/H_{2}=1$. For both
the transverse or longitudinal component, the flatness is larger for the cases
with smaller bubbles, reflecting the same trend as the PDFs of the velocity
increments at the single scale $r_{2}=10\Delta$ (see figures 12a,b or figures
12c,d). The qualitative effect of the surfactant concentration on the extreme
events as quantified by the PDF results for
$\Delta\boldsymbol{u}(r_{2}=10\Delta)$ (figures 12) can be extended to almost
all the scales. While the flatness decreases in the sequence SmTap, SmPen to
SmPen+ across the scales of the flow, the values of the three cases with
larger bubbles are very similar.

(a)

(b)

Figure 13: Normalized fourth-order transverse (a) and longitudinal (b)
structure functions, corresponding to the flatness of the velocity increments
along the horizontal direction. The horizontal lines in (a,b) indicate to the
Gaussian value of $3$ for the flatness.

##  6 Conclusions

In this work we have used experiments to investigate the effect of surfactants
on BIT using recently developed PSV and bubble identification techniques for
bubble/particle-laden flows. The experiments consider two bubble sizes, 3 mm
and 4 mm, and three different surfactant (1-Pentanol) concentrations for each
bubble size. The addition of surfactants changes the bubble shape, as well the
interface boundary condition from a quasi-free-slip to a no-slip condition.

To provide some reference cases, we first investigated how 1-Pentanol
influences single bubbles in terms of the influence on the size of the bubbles
formed in the system, the bubble aspect ratios, the bubble rise velocity and
the nature of their trajectories. We find that for the two bubble sizes
considered, an approximately saturated contamination state is reached at
$C_{\infty}\approx 1000$ ppm. The effect of the surfactant on the bubble
trajectories is similar to that observed in Tagawa et al. (2014), with helical
trajectories at lower 1-Pentanol concentrations and zigzag rising trajectories
for higher concentration for smaller bubble cases. The effect of the
surfactant was similar for the larger bubbles, although in the absence of the
surfactants, the trajectories of the larger bubbles appear chaotic.

Our results for bubble swarms show that for BIT-dominated flows the level of
anisotropy is strong in general, and not negligible even at small scales in
the flow, consistent with previous results. However, our results reveal that
for the same bubble size the flow anisotropy can be strongly enhanced by
increasing the surfactant concentration. We also investigated extreme events
in the flow by considering the normalized probability density functions of the
velocity increments at the scale of the bubble size. For the smaller bubbles,
the PDFs become increasingly non-Gaussian with when the surfactant
concentration is decreased. This is consistent with our previous finding that
the non-Gaussianity of the PDFs of the velocity increments becomes stronger as
the Reynolds number is decreased in the bubble-laden turbulent flows
considered here, the opposite of what occurs for internal intermittency for
single-phase turbulence which increases with increasing Reynolds number.
However, the dependency of the non-Gaussianity of the velocity increments on
the surfactant concentration is much weaker for the cases with larger bubbles.
An explanation for this difference is that the larger bubble cases have much
larger volume fractions compared to those for the smaller bubbles, and the
bubble wake volume is also much larger, with much higher turbulence intensity
generated by the wakes. Due to this, the velocity difference in and outside
the wake is smaller for the larger bubble cases than the smaller bubble cases,
and the intermittency is less dependent on the surfactants concentration.

While this study has focused on the impact of surfactants on the properties of
the liquid turbulence in BIT, another important topic is to investigate how
the surfactants impact the bubble clusters and bubble dispersion. These will
be investigated in a future study.

## Acknowledgements

The authors would like to acknowledge T. Ziegenhein for many technical
supports in the experimental method over the years. T.M. benefited from
fruitful discussions on this topic with J. Fröhlich, Y. Liao, P. Shi, S.
Heitkam, X. Xu and Y. Tagawa.

## Declaration of Interests

The authors report no conflict of interest.

## Data availability

The data that support the findings of this study are available from the first
author T.M. on request.

## References

  * Ahmed et al. (2020) Ahmed, Z., Izbassarov, D., Lu, J., Tryggvason, G., Muradoglu, M. & Tammisola, O. 2020 Effects of soluble surfactant on lateral migration of a bubble in a pressure driven channel flow. Int. J. Multiphase Flow 126, 103251.
  * Atasi et al. (2022) Atasi, O., Ravisankar, M., Legendre, D. & Zenit, R. 2022 The presence of surfactants controls the stability of bubble chains in carbonated drinks. arXiv preprint arXiv:2211.02253 .
  * Bachhuber & Sanford (1974) Bachhuber, C. & Sanford, C. 1974 The rise of small bubbles in water. Journal of Applied Physics 45 (6), 2567–2569.
  * Batchelor (1967) Batchelor, G. K. 1967 An introduction to fluid dynamics. Cambridge University Press.
  * Bel Fdhila & Duineveld (1996) Bel Fdhila, R. & Duineveld, P. C. 1996 The effect of surfactant on the rise of a spherical bubble at high reynolds and peclet numbers. Physics of Fluids 8 (2), 310–321.
  * Biferale et al. (2012) Biferale, L, Perlekar, P, Sbragaglia, M & Toschi, F 2012 Convection in multiphase fluid flows using lattice boltzmann methods. Phys. Rev. Lett. 108 (10), 104502.
  * Bröder & Sommerfeld (2007) Bröder, D & Sommerfeld, M 2007 Planar shadow image velocimetry for the analysis of the hydrodynamics in bubbly flows. Measurement Science and Technology 18 (8), 2513.
  * Cano-Lozano et al. (2016) Cano-Lozano, J. C., Martinez-Bazan, C., Magnaudet, J. & Tchoufag, J. 2016 Paths and wakes of deformable nearly spheroidal rising bubbles close to the transition to path instability. Phys. Rev. Fluids 1 (5), 053604.
  * Carter & Coletti (2017) Carter, D. W & Coletti, F. 2017 Scale-to-scale anisotropy in homogeneous turbulence. J. Fluid Mech. 827, 250–284.
  * du Cluzeau et al. (2022) du Cluzeau, A, Bois, G, Leoni, N & Toutant, A 2022 Analysis and modeling of bubble-induced agitation from direct numerical simulation of homogeneous bubbly flows. Phys. Rev. Fluids 7 (4), 044604.
  * du Cluzeau et al. (2019) du Cluzeau, A, Bois, G & Toutant, A 2019 Analysis and modelling of reynolds stresses in turbulent bubbly up-flows from direct numerical simulations. J. Fluid Mech. 866, 132–168.
  * Cuenot et al. (1997) Cuenot, B, Magnaudet, J & Spennato, B 1997 The effects of slightly soluble surfactants on the flow around a spherical bubble. J. Fluid Mech. 339, 25–53.
  * Drenckhan & Saint-Jalmes (2015) Drenckhan, W. & Saint-Jalmes, A. 2015 The science of foaming. Adv. Colloid Interface Sci. 222, 228–259.
  * Durst et al. (1986) Durst, F., Schönung, B., Selanger, K. & Winter, M. 1986 Bubble-driven liquid flows. J. Fluid Mech. 170, 53–82.
  * Eftekhari et al. (2021) Eftekhari, M., Schwarzenberger, K., Heitkam, S., Javadi, A., Bashkatov, A., Ata, S. & Eckert, K. 2021 Interfacial behavior of particle-laden bubbles under asymmetric shear flow. Langmuir 37 (45), 13244–13254.
  * Ellingsen & Risso (2001) Ellingsen, K. & Risso, F. 2001 On the rise of an ellipsoidal bubble in water: oscillatory paths and liquid-induced velocity. J. Fluid Mech. 440, 235–268.
  * Frisch (1995) Frisch, U. 1995 Turbulence: the legacy of AN Kolmogorov. Cambridge university press.
  * Frumkin & Levich (1947) Frumkin, A & Levich, V. G. 1947 On surfactants and interfacial motion. Zh. Fiz. Khim. 21, 1183–1204.
  * Hayashi & Tomiyama (2018) Hayashi, Kosuke & Tomiyama, Akio 2018 Effects of surfactant on lift coefficients of bubbles in linear shear flows. Int. J. Multiphase Flow 99, 86–93.
  * Hessenkemper et al. (2022) Hessenkemper, H., Starke, S., Atassi, Y., Ziegenhein, T. & Lucas, D. 2022 Bubble identification from images with machine learning methods. Int. J. Multiphase Flow 155, 104169\.
  * Hessenkemper & Ziegenhein (2018) Hessenkemper, H & Ziegenhein, T 2018 Particle shadow velocimetry (PSV) in bubbly flows. Int. J. Multiphase Flow 106, 268–279.
  * Hessenkemper et al. (2021a) Hessenkemper, H, Ziegenhein, T, Lucas, D & Tomiyama, A 2021a Influence of surfactant contaminations on the lift force of ellipsoidal bubbles in water. Int. J. Multiphase Flow 145, 103833.
  * Hessenkemper et al. (2021b) Hessenkemper, H., Ziegenhein, T., Rzehak, R., Lucas, D. & Tomiyama, A. 2021b Lift force coefficient of ellipsoidal single bubbles in water. Int. J. Multiphase Flow 138, 103587.
  * Innocenti et al. (2021) Innocenti, A., Jaccod, A., Popinet, S. & Chibbaro, S. 2021 Direct numerical simulation of bubble-induced turbulence. J. Fluid Mech. 918.
  * Lai & Socolofsky (2019) Lai, C. C. K. & Socolofsky, S. A. 2019 The turbulent kinetic energy budget in a bubble plume. J. Fluid Mech. 865, 993–1041.
  * Lance & Bataille (1991) Lance, M. & Bataille, J. 1991 Turbulence in the liquid phase of a uniform bubbly air-water flow. J. Fluid Mech. 222, 95–118.
  * Langevin (2014) Langevin, D. 2014 Rheology of adsorbed surfactant monolayers at fluid surfaces. Annu. Rev. Fluid Mech. 46, 47–65.
  * Legendre et al. (2009) Legendre, D., Lauga, E. & Magnaudet, J. 2009 Influence of slip on the dynamics of two-dimensional wakes. J. Fluid Mech. 633, 437–447.
  * Levich (1962) Levich, V. G. 1962 Physicochemical hydrodynamics. Prentice-Hall Inc.
  * Lohse (2018) Lohse, D. 2018 Bubble puzzles: From fundamentals to applications. Phys. Rev. Fluids 3 (11), 110504.
  * Lohse (2022) Lohse, D. 2022 Fundamental fluid dynamics challenges in inkjet printing. Annu. Rev. Fluid Mech. 54, 349–382.
  * Loubière & Hébrard (2004) Loubière, Karine & Hébrard, Gilles 2004 Influence of liquid surface tension (surfactants) on bubble formation at rigid and flexible orifices. Chem. Eng. Process 43 (11), 1361–1369.
  * Lu et al. (2017) Lu, J., Muradoglu, M. & Tryggvason, G. 2017 Effect of insoluble surfactant on turbulent bubbly flows in vertical channels. Int. J. Multiphase Flow 95, 135–143.
  * Lu & Tryggvason (2013) Lu, J. & Tryggvason, G. 2013 Dynamics of nearly spherical bubbles in a turbulent channel upflow. J. Fluid Mech. 732, 166–189.
  * Ma et al. (2022) Ma, T., Hessenkemper, H., Lucas, D. & Bragg, A. D. 2022 An experimental study on the multiscale properties of turbulence in bubble-laden flows. J. Fluid Mech. 936, A42.
  * Ma et al. (2020) Ma, T., Lucas, D., Jakirlić, S. & Fröhlich, J. 2020 Progress in the second-moment closure for bubbly flow based on direct numerical simulation data. J. Fluid Mech. 883, A9.
  * Ma et al. (2021) Ma, T., Ott, B., Fröhlich, J. & Bragg, A. D 2021 Scale-dependent anisotropy, energy transfer and intermittency in bubble-laden turbulent flows. J. Fluid Mech. 927, A16.
  * Maeda et al. (2021) Maeda, K., Date, M., Sugiyama, K., Takagi, S. & Matsumoto, Y. 2021 Viscid–inviscid interactions of pairwise bubbles in a turbulent channel flow and their implications for bubble clustering. J. Fluid Mech. 919, A30.
  * Magnaudet & Eames (2000) Magnaudet, J. & Eames, I. 2000 The motion of high-reynolds-number bubbles in inhomogeneous flows. Annu. Rev. Fluid Mech. 32 (1), 659–708.
  * Manikantan & Squires (2020) Manikantan, H. & Squires, T. M 2020 Surfactant dynamics: hidden variables controlling fluid flows. J. Fluid Mech. 892, P1.
  * Masuk et al. (2021) Masuk, A. U. M., Salibindla, A. K. R. & Ni, R. 2021 Simultaneous measurements of deforming hinze-scale bubbles with surrounding turbulence. J. Fluid Mech. 910, A21.
  * Mathai et al. (2020) Mathai, V., Lohse, D. & Sun, C. 2020 Bubbly and buoyant particle–laden turbulent flows. Annu. Rev. Condens. Matter Phys. 11 11, 529–559.
  * Mclaughlin (1996) Mclaughlin, John B 1996 Numerical simulation of bubble motion in water. Journal of colloid and interface science 184 (2), 614–625.
  * Mendez-Diaz et al. (2013) Mendez-Diaz, S., Serrano-García, J. C., Zenit, R. & Hernández-Cordero, J. A. 2013 Power spectral distributions of pseudo-turbulent bubbly flows. Phys. Fluids 25 (4).
  * Moore (1963) Moore, D. W. 1963 The boundary layer on a spherical gas bubble. J. Fluid Mech. 16 (2), 161–176.
  * Mougin & Magnaudet (2001) Mougin, G. & Magnaudet, J. 2001 Path instability of a rising bubble. Phys. Rev. Lett. 88 (1), 014502.
  * Mougin & Magnaudet (2006) Mougin, G. & Magnaudet, J. 2006 Wake-induced forces and torques on a zigzagging/spiralling bubble. J. Fluid Mech. 567, 185–194.
  * Néel & Deike (2021) Néel, B & Deike, L. 2021 Collective bursting of free-surface bubbles, and the role of surface contamination. J. Fluid Mech. 917, A46.
  * Palaparthi et al. (2006) Palaparthi, R., Papageorgiou, D. T. & Maldarelli, C. 2006 Theory and experiments on the stagnant cap regime in the motion of spherical surfactant-laden bubbles. J. Fluid Mech. 559, 1–44.
  * Pandey et al. (2020) Pandey, V., Ramadugu, R. & Perlekar, P. 2020 Liquid velocity fluctuations and energy spectra in three-dimensional buoyancy-driven bubbly flows. J. Fluid Mech. 884, R6.
  * Pesci et al. (2018) Pesci, C., Weiner, A., Marschall, H. & Bothe, D. 2018 Computational analysis of single rising bubbles influenced by soluble surfactant. J. Fluid Mech. 856, 709–763.
  * Rensen et al. (2005) Rensen, J., Luther, S. & Lohse, D. 2005 The effect of bubbles on developed turbulence. J. Fluid Mech. 538, 153–187.
  * Riboux et al. (2013) Riboux, G., Legendre, D. & Risso, F. 2013 A model of bubble-induced turbulence based on large-scale wake interactions. J. Fluid Mech. 719, 362–387.
  * Riboux et al. (2010) Riboux, G., Risso, F. & Legendre, D. 2010 Experimental characterization of the agitation generated by bubbles rising at high Reynolds number. J. Fluid Mech. 643, 509–539.
  * Risso (2018) Risso, F. 2018 Agitation, mixing, and transfers induced by bubbles. Annu. Rev. Fluid Mech. 50, 25–48.
  * Roghair et al. (2011) Roghair, I., Mercado, J. M., Van Sint Annaland, M., Kuipers, H., Sun, C. & Lohse, D. 2011 Energy spectra and bubble velocity distributions in pseudo-turbulence: Numerical simulations vs. experiments. Int. J. Multiphase Flow 37, 1093–1098.
  * Santarelli et al. (2016) Santarelli, C., Roussel, J. & Fröhlich, J. 2016 Budget analysis of the turbulent kinetic energy for bubbly flow in a vertical channel. Chem. Eng. Sci. 141, 46–62.
  * Schlüter et al. (2021) Schlüter, M., Herres-Pawlis, S., Nieken, U., Tuttlies, U. & Bothe, D. 2021 Small-scale phenomena in reactive bubbly flows: Experiments, numerical modeling, and applications. Annu. Rev. Chem. Biomol. Eng. 12, 625–643.
  * Shew & Pinton (2006) Shew, W. L. & Pinton, J.-F. 2006 Dynamical model of bubble path instability. Phys. Rev. Lett. 97 (14), 144508\.
  * Soligo et al. (2019) Soligo, G., Roccon, A. & Soldati, A. 2019 Breakage, coalescence and size distribution of surfactant-laden droplets in turbulent flow. J. Fluid Mech. 881, 244–282.
  * Sreenivasan & Antonia (1997) Sreenivasan, K. R. & Antonia, R. A. 1997 The phenomenology of small-scale turbulence. Annu. Rev. Fluid Mech. 29 (1), 435–472.
  * Stone (1994) Stone, H. A. 1994 Dynamics of drop deformation and breakup in viscous fluids. Annu. Rev. Fluid Mech. 26 (1), 65–102.
  * Tagawa et al. (2014) Tagawa, Y., Takagi, S. & Matsumoto, Y. 2014 Surfactant effect on path instability of a rising bubble. J. Fluid Mech. 738, 124.
  * Takagi & Matsumoto (2011) Takagi, S. & Matsumoto, Y. 2011 Surfactant effects on bubble motion and bubbly flows. Annu. Rev. Fluid Mech. 43, 615–636.
  * Takagi et al. (2009) Takagi, S., Ogasawara, T., Fukuta, M. & Matsumoto, Y. 2009 Surfactant effect on the bubble motions and bubbly flow structures in a vertical channel. Fluid dynamics research 41 (6), 065003.
  * Takagi et al. (2008) Takagi, S., Ogasawara, T. & Matsumoto, Y. 2008 The effects of surfactant on the multiscale structure of bubbly flows. Phil. Trans. R. Soc. A 366 (1873), 2117–2129.
  * Tomiyama et al. (2002) Tomiyama, A., Celata, G.P., Hosokawa, S. & Yoshida, S. 2002 Terminal velocity of single bubbles in surface tension force dominant regime. Int. J. Multiphase Flow 28 (9), 1497–1519.
  * Townsend (1949) Townsend, A. 1949 The fully developed wake of a circular cylinder. Aust. J. Chem. 2 (4), 451–468.
  * Veldhuis et al. (2008) Veldhuis, C., Biesheuvel, A. & Van Wijngaarden, L. 2008 Shape oscillations on bubbles rising in clean and in tap water. Phys. Fluids 20 (4), 040705.
  * Verschoof et al. (2016) Verschoof, R. A., Van Der Veen, R. CA, Sun, C. & Lohse, D. 2016 Bubble drag reduction requires large bubbles. Phys. Rev. Lett. 117 (10), 104502.
  * Ybert & di Meglio (1998) Ybert, C. & di Meglio, J.-M. 1998 Ascending air bubbles in protein solutions. The European Physical Journal B-Condensed Matter and Complex Systems 4 (3), 313–319.
  * Ybert & di Meglio (2000) Ybert, C. & di Meglio, J.-M. 2000 Ascending air bubbles in solutions of surface-active molecules: influence of desorption kinetics. The European Physical Journal E 3 (2), 143–148.
  * Ziegenhein & Lucas (2017) Ziegenhein, T. & Lucas, D. 2017 Observations on bubble shapes in bubble columns under different flow conditions. Experimental Thermal and Fluid Science 85, 248–256.
  * Ziegenhein & Lucas (2019) Ziegenhein, T. & Lucas, D. 2019 The critical bubble diameter of the lift force in technical and environmental, buoyancy-driven bubbly flows. Int. J. Multiphase Flow 116, 26–38.