1. Introduction
In fully developed turbulent flows, universality is observed for velocity and temperature fluctuations for the second-order structure functions and the corresponding spectra. According to the Kolmogorov–Obukhov–Corrsin theory (Kolmogorov Reference Kolmogorov1941; Obukhov Reference Obukhov1949; Corrsin Reference Corrsin1951) and verified extensively in experiment and numerical simulations, the universal scaling of the spectrum of a scalar exhibits a 
$-5/3$ scaling (
$E_{\theta }(k) \propto k^{-5/3}$) in the inertial-convective range for high enough Reynolds and Péclet number (see e.g. Monin & Yaglom Reference Monin and Yaglom1975). Here 
$k$ is the wavenumber and 
$\theta$ is a passive scalar (Monin & Yaglom Reference Monin and Yaglom1975) such as the concentration of a chemical substance or the temperature, in case that buoyancy is not relevant. If additional complexity is added to the system, different scaling phenomena can appear, such as having an active scalar introducing buoyancy effects (e.g. Rayleigh–Bénard convection) (Lohse & Xia Reference Lohse and Xia2010), adding viscoelasticity to the carrier liquid (Steinberg Reference Steinberg2020), or by adding finite particles (Risso Reference Risso2018; Brandt & Coletti Reference Brandt and Coletti2022) or bubbles (Balachandar & Eaton Reference Balachandar and Eaton2010; Lohse Reference Lohse2018; Risso Reference Risso2018; Mathai, Lohse & Sun Reference Mathai, Lohse and Sun2020). In this work, we do the latter and introduce a large number of finite-size rising bubbles to the liquid flow (Balachandar & Eaton Reference Balachandar and Eaton2010; Lohse Reference Lohse2018; Risso Reference Risso2018; Mathai et al. Reference Mathai, Lohse and Sun2020). Large rising bubbles are defined here as having a bubble Reynolds number 
${Re}_{{bub}}=V_r d /\nu = {O}(10^2)$, where 
$V_r$ is the bubble rise velocity relative to the mean vertical velocity of the carrier liquid, 
$d$ is the mean area-equivalent bubble diameter (as compared with a spherical bubble) and 
$\nu$ is the kinematic viscosity of the carrier liquid.
First observed and theoretically addressed in Lance & Bataille (Reference Lance and Bataille1991), the introduction of bubbles lead to the emergence of a 
$k^{-3}$ scaling of the velocity fluctuation (energy) spectra 
$E_{u}(k)$. They argued that in a statistical steady state, the 
$k^{-3}$ scaling in the spectrum can be achieved by balancing the viscous dissipation and the energy production due to the rising bubbles in the spectral space (Lance & Bataille Reference Lance and Bataille1991). The 
$k^{-3}$ scaling can also result from spatial and temporal velocity fluctuations in homogeneous bubbly flows (Risso Reference Risso2018), and is found to be robust and is observed for a homogeneous bubble swarm for a wide range of parameters: 
$10 \leq {Re}_{{bub}} \leq 1000$ and 
$1 \leq {We} \leq 4$, where 
${We}$ is the Weber number (Risso Reference Risso2018; Pandey, Ramadugu & Perlekar Reference Pandey, Ramadugu and Perlekar2020). Recently, Ma et al. (Reference Ma, Hessenkemper, Lucas and Bragg2022) employed particle shadow velocimetry in bubble-laden turbulent flow and obtained a high-resolution two-dimensional velocity field. They revealed longitudinal and traverse velocity structure functions up to twelfth order, in which the second-order statistics (one-dimensional longitudinal energy spectra) also show a 
$-3$ scaling. In contrast, point particle simulations do not show 
$k^{-3}$ scaling (Mazzitelli & Lohse Reference Mazzitelli and Lohse2009), because the wakes behind the bubbles which are obviously absent for point particles are crucial for the emergence of a 
$-3$ scaling (in either frequency or wavenumber space).
Although there are abundant studies on the energy spectra in high-
${Re}$ bubbly flow, the study of scalar spectra in such flow is limited. Alméras et al. (Reference Alméras, Cazin, Roig, Risso, Augier and Plais2016) investigated the time-resolved concentration fluctuations of a fluorescent dye (passive scalar) in a confined bubbly thin cell in which the scalar spectral scaling 
$f^{-3}$ is observed, where 
$f$ is the frequency. For the three-dimensional case, Gvozdić et al. (Reference Gvozdić, Alméras, Mathai, Zhu, van Gils, Verzicco, Huisman, Sun and Lohse2018) investigated the scalar spectra in a bubbly column in a vertical convection set-up, but they were unable to resolve the frequencies related to the 
$-3$ scaling of the energy spectra in bubble-induced turbulence (Gvozdić et al. Reference Gvozdić, Alméras, Mathai, Zhu, van Gils, Verzicco, Huisman, Sun and Lohse2018). The obvious downside of using dyes is that the water gets contaminated and only very short runs can be performed before the water has to be replaced completely. By using temperature as the approximately passive scalar, and the active cooling of the set-up during the return stage, it is possible to maintain an overall constant temperature of the set-up, allowing for continuous measurements for hours or days if needed.
In the present experimental work, by utilizing a fast-response temperature probe, we are able to capture the gradual change of the passive scalar (thermal) spectral scaling for increasing bubble concentration 
$\alpha$ in a bubbly turbulent thermal mixing layer. In addition, we investigate the transition scale of the 
$-3$ scaling as it illuminates the physical origin of the 
$-3$ subrange of both the scalar and the energy spectra. The transition frequency to the 
$-3$ subrange for energy spectra indeed has yet to be exactly identified (Risso Reference Risso2018).
2. Experimental set-up and methods
2.1. Experimental facility
We utilise the Twente Mass and Heat Transfer Tunnel (Gvozdić et al. Reference Gvozdić, Dung, van Gils, Bruggert, Alméras, Sun, Lohse and Huisman2019) which creates a turbulent vertical channel flow, and allows bubbles to be injected and the liquid to be heated in a controlled manner. The turbulence is actively stirred using an active grid, see figure 1(a). Using this facility we create a turbulent thermal mixing layer (de Bruyn Kops & Riley Reference de Bruyn Kops and Riley2000) in water with a mean flow velocity of approximately 0.5 m s
$^{-1}$, resulting in a Reynolds number of 
${Re} = 2 \times 10^4$. The global gas volume fraction 
$\alpha$ in the measurement section is measured by a differential pressure transducer of which the two ends are connected to the top and bottom of the measurement section. A fast-response thermistor (Amphenol Advanced Sensors type FP07, with a response time of 7 ms in water) is placed into the middle of the measurement section in order to measure the temperature. An AC bridge with a sinusoidal frequency of 1.3 kHz drives the thermistor in order to reduce electric noise. The bridge potential is measured by a lock-in amplifier (Zurich Instruments MFLI) with the acquisition frequency of 13.4 kHz and such that the cut-off frequency that is used is (13.4/2) kHz, which is far from the effective range of the thermistor which is limited by the thermal capacity of the probe. A temperature-controlled refrigerated circulator with water bath (PolyScience PD15R-30) with a temperature stability of 
$\pm$0.005 K is used for calibrating the thermistor. We employ the temperature-resistance characteristic equation proposed by Steinhart & Hart (Reference Steinhart and Hart1968) over a short range of temperature of interest (in which the terms are kept only up to the first order to avoid over-fitting) in order to convert the measured resistance of the thermistor to the corresponding temperature. The thermistor has a relative precision of 1 mK. With relative we mean here relative to a baseline temperature, i.e. a temperature difference (
${\rm \Delta} T$), as we are mainly interested in temperature changes, which have a different precision than 
$T$ itself. The typical large-scale temperature difference from the heated side to the non-heated side at the mid-height of the measurement section of height 
$h = 1$ m is about 0.2 K. The working liquid is decalcified tap water at around 
$22.7\,^\circ$C, with a Prandtl number 
${Pr}=\nu /\kappa =6.5$, where 
$\nu$ and 
$\kappa$ are the kinematic viscosity and thermal diffusivity of water, respectively.
(a) Schematic of the experimental set-up used in the current study. Bubbles are injected by 120 needles (green). The set-up is equipped with 12 heaters. A power of 2250 W is supplied to the left half of the heaters (red), whereas the right half (white) is left unpowered, creating a turbulent thermal mixing layer. Each heating cartridge (Watlow, Firerod J5F-15004) has a diameter of 12.7 mm. The centre-to-centre distance between two consecutive heaters is 37.5 mm. The turbulence is stirred by an active grid with a total of 15 engines connected to diamond flaps (yellow) (Gvozdić et al. Reference Gvozdić, Dung, van Gils, Bruggert, Alméras, Sun, Lohse and Huisman2019). The measurement section (blue area) has a height of 1.00 m, a width of 0.3 m, and is 0.04 m thick. (b–d) High-speed snapshots of the bubble flow, as seen from the side for the volume fractions 
$\alpha = 0.3\,\%$, 
$1.4\,\%$ and 
$4.7\,\%$. Scale bars of 20 mm are added to each figure. (e) Distribution charts of the bubble diameter obtained using image analysis of high-speed footage for a variety of gas volume fractions 
$\alpha$. The resolution of the footage prohibits us from measuring bubbles (dashed line) of size 
$<1.4~mm$. Dots indicate the mean bubble diameter 
$d$. (f) Histograms of the bubble aspect ratio 
$\chi$ obtained by fitting the bubbles on the high-speed footage with ellipses. Here 
$\chi$ is defined as the ratio of the semi-major axis divided by the semi-minor axis. Each black dot represents the mean aspect ratio 
$\chi$ and each error bar represent one standard deviation of the distribution. See also table 2 for the details of the bubble properties.
For the velocity measurements, we employ laser Doppler anemometry (LDA) (Dantec) and constant temperature anemometry (CTA) using a hot film sensor (Dantec 55R11, 70 
$\mathrm {\mu }$m diameter with a 2 
$\mathrm {\mu }$m nickel coating, 1.3 mm long and up to 30 kHz frequency response) placed in the middle of the tunnel. CTA is used for calculating the velocity spectra and the corresponding method employed in this work can be found in Appendix A. We note that because a CTA measurement requires negligible temperature variation of the flow, the velocity is measured when the heaters are switched off, but left in place such as to not alter the incoming flow. The velocity statistics are not influenced by the temperature because the fluctuations in the flow are completely dominated by the turbulence generated by the active grid and the upward mean flow rather than by buoyancy caused by density differences due to the temperature differences, i.e. it is an approximately passive scalar. To substantiate this statement, we must compare the relative importance of the buoyancy force (parallel to the streamwise direction) to the inertial force due to the advection of a passive scalar. The ratio of these forces equals the ratio of the Grashof number 
${Gr}$ and the square of the Reynolds number 
${Re}_{0.5h}^2$ (Schlichting & Gersten Reference Schlichting and Gersten2017). Here we base these numbers on the vertical distance between 
$\pm 0.25h$ from the middle of the measurement section, the temperature difference across such height 
$\varDelta _{0.5h} \leq 50$ mK and the mean liquid velocity 
$U = 0.5$ ms
$^{-1}$. This gives 
${Gr}/{Re}_{0.5h}^2 =0.5hg\beta \varDelta _{0.5h}/U^2 \leq {O}(10^{-3})$, where 
$g$ and 
$\beta$ are the gravitational constant and volumetric expansion coefficient of water, respectively. The small value of this ratio implies that the temperature can be seen as passive scalar under our flow conditions.
2.2. Measurement conditions: single-phase turbulent characterisation
Next, we present the turbulent characteristics in the single-phase turbulent thermal mixing layer. First, we characterise the velocity field. Here, LDA (Dantec) is used for measuring the horizontal (
$u_x$) and vertical (
$u_z$) velocity in order to characterise the velocity fluctuations. Throughout the paper we use a prime (
$'$) to denote the standard deviation of any quantity, and therefore the standard deviation of the velocity in the single-phase case is 
$u_0'$ and we approximate 
$u_0'\approx u_z'$. In the single-phase case, at the middle of the tunnel, 
$u_x'/u_z' = 0.69 \pm 0.07$, where the statistical uncertainty is one standard deviation. This is not ideally isotropic. Therefore, for single-phase case, the estimation of the relevant scales using the assumption of isotropy in the following is only a rough approximation. The viscous dissipation rate 
$\epsilon$ is estimated by finding the longitudinal second-order structure function 
$D_{LL}(r)$ in single-phase flow, where 
$r$ is the spatial distance between two points (measured by CTA and used Taylor's frozen flow hypothesis) and using the velocity two-third law 
$D_{LL}= C_2 (\epsilon r)^{2/3}$ which is valid in the inertial range, assuming isotropy and where 
$C_2 = 2.0$ is a universal constant (Pope Reference Pope2000). As the inertial range has limited extension (see figure 2a), we first identify the value of 
$r$ such that 
$D_{LL}(r)/r^{2/3}$ becomes locally flat versus 
$r$, and then employ the velocity two-third law. We then obtain the Kolmogorov length scale 
$\eta = (\nu ^3/\epsilon )^{1/4} = 0.22$ mm and the dissipation time scale 
$\tau _{\eta } = (\nu /\epsilon )^{1/2} = 49$ ms. Moreover, we estimate the Taylor-microscale 
$\lambda _u = \sqrt {15 \nu (2k/3)/\epsilon } = 4.9$ mm and the Taylor–Reynolds number 
$Re_{\lambda _u}=(2k/3)^{1/2} \lambda _u/\nu = 130$, where 
$k=3u_0'/2$ is the average turbulent kinetic energy, by assuming isotropy. To estimate the velocity integral length scale 
$L_u$, we use 
$\epsilon = C_{\epsilon } u_0'^3/L_u$, where 
$C_{\epsilon }=0.9$ (Valente & Vassilicos Reference Valente and Vassilicos2012; Vassilicos Reference Vassilicos2015).
The power spectra 
$P_i(f)$ normalised by their respective variance of (a) the velocity fluctuations (
$P_u(f)$) and of (b) the temperature fluctuations (
$P_{\theta }(f)$) for 
$\alpha$ from 0 % (purple) to 5.2 % (red). The compensated spectra 
$f^3P_{i}(f)$ and 
$f^{5/3}P_{\theta }(f)$ are shown as insets. The scalar spectra are cut at the limit of 143 Hz, corresponding to the highest resolved frequency of the thermistor.
Second, for the scalar field characterisation, the temperature fluctuations are measured by a thermistor. It is well-known that passive scalar fields are not isotropic at small scales when there is a large-scale mean temperature gradient (Warhaft Reference Warhaft2000) but for the order of estimation of the quantities in the following, we assume isotropy. First, to estimate the scalar dissipation rate 
$\epsilon _{\theta }$. To do so, we again assume isotropy and use the temperature (scalar) two-third law 
$D_{\theta \theta }(r) = \widetilde {C_{\theta }}\epsilon _{\theta }\epsilon ^{-1/3} r^{2/3}$ which is valid in the inertial-convective subrange (Monin & Yaglom Reference Monin and Yaglom1975), where 
$D_{\theta \theta }$ is the second-order structure function of the temperature (scalar) fluctuation and 
$\widetilde {C_{\theta }} \approx 0.25^{-1} C_{\theta }$ (Monin & Yaglom Reference Monin and Yaglom1975) with 
$C_{\theta } \sim 0.5$ the Obukhov–Corrsin constant (Mydlarski & Warhaft Reference Mydlarski and Warhaft1998; Warhaft Reference Warhaft2000). Again we identify the plateau of the plot 
$D_{\theta \theta }(r)/r^{2/3}$ versus 
$r$ in order to locate the position of 
$r$ at which we employ the temperature two-third law. Second, we estimate the scalar Taylor microscale 
$\lambda _{\theta } = \sqrt {6 \kappa \langle T'^2 \rangle /\epsilon _{\theta }} =2.1$ mm (see e.g. Yasuda et al. Reference Yasuda, Gotoh, Watanabe and Saito2020), again assuming isotropy. We find that our turbulent thermal mixing layer, for the single-phase case, has Péclet numbers (based on both the velocity and temperature Taylor microscales) 
${Pe}_{\lambda _u} = u_0'\lambda _u/\kappa = 870$ and 
${Pe}_{\lambda _{\theta }}=u_0'\lambda _{\theta }/\kappa = 520$, respectively (Yasuda et al. Reference Yasuda, Gotoh, Watanabe and Saito2020). The traditional choice of the microscale is 
$\lambda _u$ (Mydlarski & Warhaft Reference Mydlarski and Warhaft1998) but numerical simulations have uncovered that the small-scale anisotropy of the scalar field depends on the Péclet number based on 
$\lambda _{\theta }$ (Yasuda et al. Reference Yasuda, Gotoh, Watanabe and Saito2020). Thus, we give both values here. For the length scale of the scalar dissipation, because 
${Pr} = 6.5 > 1$, we use the Batchelor scale to characterise the smallest scale of the passive scalar field which is given by 
$\eta _{\theta } = \eta {Pr}^{-1/2}$ (Tennekes & Lumley Reference Tennekes and Lumley1972). As the mean temperature profile for the single-phase turbulent mixing layer follows a self-similar error-function profile (Ma & Warhaft Reference Ma and Warhaft1986; de Bruyn Kops & Riley Reference de Bruyn Kops and Riley2000), we fitted the profile accordingly and obtain the large-scale characteristic length 
$L_{\theta }$ of the single-phase thermal mixing layer from the fitting parameters (see Appendix B for details). The turbulent characteristics of the single-phase turbulent mixing layer are summarised in table 1.
Relevant turbulent flow parameters of the single-phase turbulent thermal mixing layer in the current study (
$\alpha = 0\,\%$, mean liquid flow speed 
$= 0.5$ m s
$^{-1}$). Note that isotropy is assumed for the estimation of the above quantities except 
$L_{\theta }$. For the definitions of the parameters and the discussion on the isotropy, see § 2.
2.3. Measurement conditions: bubble properties
Bubbles are created by injecting regular air through needles, see figure 1(a). Using high-speed imaging in a backlit configuration we characterise the size of our bubbles for a variety of volume fractions, see the example photos in figure 1(b–d). We analyse the data by making use of the circular Hough transform to detect the individual (potentially overlapping) bubbles. We track the bubbles over time such as to obtain the mean bubble rise velocity, see table 2. In figure 1(e) we show the probability density function (PDF) of the bubble diameter for several 
$\alpha$ in the form of distribution charts. The mean bubble diameter obtained from the image analysis is approximately constant but there are changes in the bubble distribution over 
$\alpha$.
Experimental parameters as a function of 
$\alpha$: velocity fluctuations of the liquid 
$U' \equiv \sqrt {2 u_x'^2 + u_z'^2 }$ with primed velocities referring to the standard deviation of these velocity components, area-equivalent bubble diameter 
$d$, the Weber number based on the turbulent velocity 
$U'$ for the bubbles 
${We}_{U'}=\rho U'^2 d/\gamma$ (with 
$\gamma$ the surface tension), the aspect ratio of the bubble diameter 
$\chi$, the bubble rise velocity 
$V_{bub}$, the bubble relative (to the liquid phase) rise velocity 
$V_r \equiv V_{{bub}} - \langle u_z \rangle$, the mean Weber number based on the bubble relative rise velocity 
${We}_{V_r}=\rho V_r^2 d/\gamma$, the bubble Reynolds number 
${Re}_{bub} = V_r d/\nu$, the fitting parameters 
$f_t$ and 
$f_L$ (below (4.1)) and 
$f_{c,A}$ (4.3) (Alméras et al. Reference Alméras, Mathai, Lohse and Sun2017). The bubble sizes and velocity are obtained by image analysis (see § 2). *Note that for 
$\alpha = 0.3\,\%$, we slightly overestimate the mean bubble diameter and underestimate the width of the bubble size distribution because of the resolution limit of the image analysis for the bubble diameter. The aspect ratio of the bubbles are obtained by fitting the bubbles with ellipses. The uncertainty of 
$\alpha$ stems from the precision of the differential pressure gauge (Gvozdić et al. Reference Gvozdić, Dung, van Gils, Bruggert, Alméras, Sun, Lohse and Huisman2019). For 
$U'$, 
$d$, 
${We}_{U'}$, 
$\chi$ and 
$V_{bub}$, we show the mean values and the corresponding standard deviations for the distributions. For 
$V_r$, 
${We}_{V_r}$ and 
${Re}_{bub}$, we show the corresponding mean values and the standard errors of the mean. For 
$f_L$ and 
$f_t$, we show the mean values and the corresponding confidence intervals for a 95 % confidence level obtained from the fitting. For 
$f_{c,A}$, we show the mean values and the corresponding confidence intervals for a 95 % confidence level.
First, we discuss the change of the bubble distribution for different 
$\alpha$ and the possible contributions that lead to such change. For 
$\alpha = 0.3\,\%$, when measuring the bubble diameter through image analysis, there are unresolved bubbles with sizes smaller than 1.4 mm, thus overestimating the mean bubble diameter for that small 
$\alpha = 0.3\,\%$. The larger spread of the bubble size distribution for 
$\alpha =0.3\,\%$ may be accounted for as follows. In order to produce bubbles in our set-up, the array of needles at the bottom of the set-up is fed by pressurised air from only one side (see figure 1(a) and also Gvozdić et al. Reference Gvozdić, Dung, van Gils, Bruggert, Alméras, Sun, Lohse and Huisman2019). The case for 
$\alpha = 0.3\,\%$ is close to the minimum achievable gas volume fraction at which the pressure that pushes the air out from the array of needles may become more uneven at different needle positions (the needles that are closer to the inlet valve have larger pressure) than for the higher 
$\alpha$ cases. Thus, the bubbles come out with more diverse sizes than for the cases with higher 
$\alpha$. Moreover, higher 
$\alpha$ also means that there are more bubble–bubble interactions which include more coalescence of bubbles, thus contributing to the smaller spread of the bubble sizes than what happens at lower 
$\alpha$ cases. For 
$\alpha >0.3\,\%$, the mean bubble diameter is essentially unchanged, having the value of 
$d\approx 2.5$ mm, see also table 2, column 
$d$. However, even for 
$\alpha > 0.3\,\%$, the standard deviation of the bubble sizes changes slightly, in particular dropping from 0.6 to 0.4 mm when 
$\alpha$ increases from 
$2.2\,\%$ to 
$3.0\,\%$. The distributions remain to have similar shape for higher 
$\alpha$ except having a slight increase of width at the highest 
$\alpha = 4.7\,\%$ (table 2). However, this slight change of the width of the distribution for 
$\alpha > 0.3\,\%$ can be considered minor if we see the width of the bubble distribution is around 0.5 mm for 
$\alpha > 0.3\,\%$.
For the slight decrease of the width of the bubble size distribution when 
$\alpha$ increases further from 
$0.3\,\%$, we can attribute to (1) the possibility that more even pressure distribution over the different needles that yields more evenly distributed sizes of bubbles when the bubbles comes out, and (2) more bubble–bubble interactions and coalescence of bubbles, which makes the spread of the sizes lower. There are also (3) the effect of incidence turbulence and (4) the effect of the cutting of the bubbles due to active grid (discussed in the next paragraph) which can also affect the bubble size distribution.
Second, we compare the variation of the mean bubble diameter in our case with other studies and discuss a possibility that leads to our observation. Although the mean bubble diameter remains roughly constant over 
$\alpha$, other research, which also makes use of capillary needles to inject bubbles, indicates that the bubble diameter increases with 
$\alpha$. This includes a homogeneous bubbly flow in quiescent liquid (Colombet et al. Reference Colombet, Legendre, Risso, Cockx and Guiraud2015) and a rising bubbly flow with a background upward liquid flow which is weakly turbulent (Roig & Larue de Tournemine Reference Roig and Larue de Tournemine2007). Bubble swarms with a background upward liquid flow that is strongly turbulent were investigated in Alméras et al. (Reference Alméras, Mathai, Lohse and Sun2017) and Prakash et al. (Reference Prakash, Mercado, van Wijngaarden, Mancilla, Tagawa, Lohse and Sun2016) by using the Twente Water Tunnel (Poorte & Biesheuvel Reference Poorte and Biesheuvel2002), and they find also a small increase of 
$d$ with 
$\alpha$ but the investigated range of 
$\alpha$ was limited (
$\alpha < 1\,\%$). It is not clear why the bubble sizes are roughly constant in our case. The active turbulent grid in our set-up has a similar design but with different aspect ratio and has a smaller overall size than that used in the Twente Water Tunnel (Prakash et al. Reference Prakash, Mercado, van Wijngaarden, Mancilla, Tagawa, Lohse and Sun2016; Alméras et al. Reference Alméras, Mathai, Lohse and Sun2017). It is also reported in Alméras et al. (Reference Alméras, Mathai, Lohse and Sun2017) and Prakash et al. (Reference Prakash, Mercado, van Wijngaarden, Mancilla, Tagawa, Lohse and Sun2016) that the fragmentation of bubbles caused by the active grid leads to smaller bubble sizes for high liquid mean flow. In our case, although there is a possibility that bubble diameter may increase with 
$\alpha$ before the bubbles entering the active turbulent grid, possibly the rotating speed of the active turbulent grid is too high that it cuts the bigger bubbles into smaller bubbles with sizes similar to those at lower 
$\alpha$.
Third, we obtained the aspect ratio 
$\chi$ of the bubbles by fitting the bubbles with ellipses in the high-speed footage. The bubble diameter analysis has a resolution limit of 1.4 mm, whereas the aspect-ratio analysis uses all sizes (though the tiny bubbles of 
${O}(1)$ pixels are still not selected for analysis). Bubbles with smaller sizes, in particular for 
$\alpha = 0.3\,\%$, are more spherical bubbles because of the stronger effect of surface tension. For each 
$\alpha$, we use 
${O}(200)$ bubbles for fitting to obtain the aspect ratios statistics. Figure 1(f) shows 
$\chi$ versus 
$\alpha$. The trend shows the mean be to almost monotonically increasing with 
$\alpha$, saturated at high 
$\alpha$. When 
$\alpha$ increases from 
$0.3\,\%$ to 
$1.0\,\%$, there is a considerable increase of aspect ratio 
$\chi$ from 1.4 to 1.7, see table 2. Increasing 
$\alpha$ further, the aspect ratio saturates. In attempt to explain the monotonic increase of 
$\chi$ versus 
$\alpha$, one may check the Weber number and bubble Reynolds number, which are the dimensionless parameters for a freely rising single bubble. Although there is a slight increase of 
$d$ for 
$\alpha$ being from 
$0.3\,\%$ to 
$0.6\,\%$, the relative rise velocity of the bubbles decreased dramatically when 
$\alpha$ increases. Therefore, the bubble Reynolds number drops dramatically (table 2). Moreover, a lower Weber number based on the relative rise velocity 
${We}_{V_r}$ is also a result (table 2). However, the Weber number based the liquid velocity fluctuations 
${We}_{U'}$ increases because the liquid velocity agitation dramatically increases with 
$\alpha$. There are also higher (1) bubble–bubble interactions as 
$\alpha$ increases and we have to take into account of (2) incident turbulence in our case and its interaction with the bubbles. Then, apart from the inhomogeneity of the pressure distribution at the needles that produces gas bubbles for very low 
$\alpha$, these two contributions may also affect the value of 
$\chi$ and subject to further systematic investigation. We note that in the case of freely rising bubble swarm simulation (see figure 5.12 in Roghair Reference Roghair2012), the bubbles with the same diameter become more spherical when 
$\alpha$ increases, as opposed to what we observed here. This simulation result may limit to the importance of other possible effects such as incidence turbulence as described above. Finally, we also point out that the change of shape is linked to 
$\alpha$ and we cannot exclude that the change of aspect ratio is more important than the change of 
$\alpha$ itself. We use these bubbles to provide additional stirring to the liquid.
For the cases with bubbles, because it is highly anisotropic between vertical and azimuthal directions, we estimate the velocity fluctuations of the liquid as 
$U' \equiv \sqrt {2\langle u_x'^2 \rangle + \langle u_z'^2 \rangle }$, where we assume 
$\langle u_x'^2 \rangle = \langle u_y'^2 \rangle$. We first explore the parameter space by increasing 
$\alpha$ from 
$0.0\,\%$ to 
$5.2\,\%$ for a fixed flow velocity of 0.5 m s
$^{-1}$ and with our active grid turned on. Finally, the bubble Reynolds number 
${Re}_{{bub}}=V_r d/\nu$ as obtained from visualisation is in the range between 480 (
$\alpha =0.6\,\%$) and 220 (
$\alpha =4.7\,\%$). These properties of the bubbles are summarised in table 2.
3. The energy and scalar spectra
As the scalar in the present study is advected passively by the liquid velocity, the corresponding energy spectra first need to be examined in order to understand the passive scalar spectra. In figure 2 we show the energy and scalar power spectra 
$P_u$ and 
$P_\theta$ for different gas concentrations 
$\alpha$ normalised by their corresponding variance, corresponding to the normalised energy and scalar power spectra, respectively. The methods of calculating the spectra of the velocity and temperature fluctuations for two-phase flows are discussed in Appendix A. The energy spectra for single-phase shows a limited inertial range. However, whenever 
$\alpha > 0$, a pronounced 
$-3$ scaling is observed with the presence of an energy ‘bump’ just before this 
$-3$ scaling when 
$\alpha$ is large enough, which is consistent with the results of Alméras et al. (Reference Alméras, Mathai, Lohse and Sun2017). We find that such 
$-3$ scaling occurs roughly between 100 and 1000 Hz and has a higher energy content than the single phase in the same frequency range. The physical mechanism of such small-scale enhanced velocity fluctuations is attributed to bubble-induced agitations, in which the velocity fluctuations produced by the bubbles are directly dissipated by viscosity (Lance & Bataille Reference Lance and Bataille1991; Prakash et al. Reference Prakash, Mercado, van Wijngaarden, Mancilla, Tagawa, Lohse and Sun2016; Risso Reference Risso2018).
Next, we examine the passive scalar (thermal) spectra. A scaling of 
$-5/3$ for 
$1\ {\rm Hz} \leq f \leq 15$ Hz is observed for the single-phase scalar spectrum. The phenomenon that a developed inertial range is observed in a passive scalar spectrum but not in the corresponding energy spectrum due to limited Reynolds number was found previously by Jayesh, Tong & Warhaft (Reference Jayesh, Tong and Warhaft1994). As opposed to the relative high-frequency (small-scale) fluctuation enhancement in the energy spectra, the scalar spectra show relative high-frequency (small-scale) diminution when 
$\alpha$ increases from zero. The 
$-5/3$ scaling is followed by a steeper slope when 
$\alpha$ increases but such change of scaling is gradual and approximately saturated to a rough scaling of 
$-3$ for high enough 
$\alpha$ (
$\alpha \approx 4.7\,\%$) around 
$f={O}(10\ \text {Hz})$. However, as seen also in the inset of figure 2(b) for 
$-5/3$ compensated spectra or the spectral local exponent in figure 3(b), the original 
$-5/3$ scaling becomes less pronounced and has a steeper slope also at the frequency range before the developed 
$-3$ scaling for the thermal spectra. The main differences between the energy and the passive scalar spectra are that the rising bubbles produce additional fluctuations in the former but diminish the fluctuations in the latter. As we explain in the following, this can be understood due to the added small-scale mixing of the thermal fluctuations by the bubbles which smoothes the temperature field. The emergence of the 
$-3$ scaling is abrupt in the former but gradual in the latter; and the 
$-3$ scaling occurs at 
${O}(100\ \text {Hz})$ for the former but 
${O}(10 \text {Hz})$ for the latter. Clearly, there are differences in the physical mechanisms that lead to the same 
$-3$ scaling for the energy and passive scalar spectra. We first clarify the physical and mathematical background before speculating on the emergence of the 
$-3$ scaling in the scalar spectra.
(a) The power spectral densities 
$P_{\theta }(f)$ (normalised by their respective variance) of the temperature fluctuation from the measurements (dots) and the corresponding fits (4.1) in logarithmic scale for the gas volume fractions 
$\alpha = 1.0\,\%$ and 
$\alpha = 5.2\,\%$ at the frequency range of 1.2–15 Hz (colours are the same as in figure 2), plotted together to the purpose of showing the quality of the fits. The two grey sloped lines labelled with ‘
$-5/3$’ and ‘
$-3$’ indicate the scaling behaviours of the two limits. The vertical grey lines labelled by 
$f_t$ indicate the fitted transition frequency for the two cases. (b) The local scaling exponent 
$\xi (f)$ (defined by (4.2)) of 
$P_{\theta }(f)$ over an interval of more than three decades for different gas volume fractions 
$\alpha$. Two grey horizontal lines indicates the slopes 
$-5/3$ and 
$-3$. The details of obtaining the local slope is described in § 4. (c) Comparison of 
$f_{c,A}(\alpha )$ (open symbols, (4.3)) and 
$f_t(\alpha )$ (solid symbols). The uncertainty of 
$\alpha$ is related to the precision of the differential pressure gauge (Gvozdić et al. Reference Gvozdić, Dung, van Gils, Bruggert, Alméras, Sun, Lohse and Huisman2019), whereas the error bars of 
$f_d$ and 
$f_{c,A}$ are the confidence intervals from the fitting (4.1) for a 
$95\,\%$ confidence level. The fitting results of 
$f_t$ are also listed in table 2.
The main commonalities and differences of the physical settings of the velocity and temperature fields are as follows. One common aspect is that along the measurement section, both the viscosity and diffusivity reduce the velocity and temperature fluctuations, respectively. On the other hand, the mean velocity field is homogeneous in the bulk of the measurement section whereas the temperature field has a large-scale mean temperature gradient. The high-
${Re}$ bubbles have viscous boundary layers at their air–water interface but the thermal boundary layer is expected to be minimal as the heat capacity of the bubbles is tiny, and their thermal response is fast with at most a time scale of (conductive case, no thermal convection) 
$t=\delta ^2/D_{air}={O}(50\ {\rm ms})$, where 
$\delta = 0.4d$ is the penetration depth for a spherical bubble and 
$D_{air}$ the thermal diffusivity of air. If convective flow inside the bubble is included this typical time-scale is drastically reduced, implying that the bubble can adjust its temperature very quickly to its surroundings. Moreover, if we compare the heat capacities of the water and the bubble we obtain 
$(95\,\% C_{p,{water}}\rho _{{water}})/(5\,\% C_{p,{air}}\rho _{{air}}) \approx 65{\,}000$, meaning that the bubbles cannot be very effective carriers of thermal energy as compared with the water. The viscous boundary layers transport momentum to the liquid phase at the scale related to bubbles, producing velocity agitations. The absence of thermal boundary layers implies no heat transport is present at the scale related to bubbles. Therefore, the only source of scalar fluctuations is from the mean temperature gradient whereas the source of velocity fluctuations is from bubble agitations.
Mathematically, for the spectral behaviour of a passive scalar, one can consider the general spectral equation for the three-dimensional scalar spectra 
$E_{\theta }(k,t)$ derived from the advection–diffusion equation (Monin & Yaglom Reference Monin and Yaglom1975) which neglects viscous heating,
\begin{equation} \frac{\partial}{\partial t} E_{\theta}(k,t) + 2\kappa k^2 E_{\theta}(k,t) = T_{\theta}(k,t) + \varPi_{\theta}(k,t), \end{equation}
where 
$t$ is the time and 
$T_{\theta }$ and 
$\varPi _{\theta }$ are the local net transfer and the production at wavenumber 
$k$, respectively. One can employ Taylor's frozen-flow hypothesis to transform from the frequency domain to the wavenumber domain. In the present case, for increasing 
$\alpha$, the single-phase 
$-5/3$ scaling of the spectrum gradually transitions to a 
$-3$ scaling at higher 
$k$ (or 
$f$).
With the previously discussed physical settings and scaling behaviours, we now speculate on the physical mechanisms on the 
$-3$ scaling of the scalar spectra in bubbly flows. The 
$-5/3$ range is considered to be the inertial range where there is negligible net local transfer, negligible production and neglected diffusivity (and viscosity) effect (Pope Reference Pope2000). The scalar production due to the mean temperature gradient occurs at a smaller frequency than the 
$-5/3$ frequency range (lower than 1 Hz). As the production due to the presence of bubbles is negligible as discussed previously, we speculate that 
$\varPi _{\theta } \approx 0$ at the scalar 
$-3$ subrange. Furthermore, we speculate that the rising bubbles enhance the mixing of temperature by homogenising the temperature with increasing 
$\alpha$ at the scale related to the bubbles and, thus, we observe a faster drop (steeper slope) in the passive scalar spectra. When 
$\alpha$ is large enough, the smoothing of the temperature inhomogeneity due to the bubbles is so strong that the thermal fluctuations are directly smoothed out by the molecular thermal diffusivity which causes the saturation of the scaling. In other words, the thermal fluctuations that originally have cascading behaviour (as in the single-phase case) passing them from larger scales to smaller scales and ceasing by the molecular thermal diffusion. When this situation happens, in this regime, we speculate that 
$T_{\theta } = T_{\theta } (\epsilon _{\theta },k)$, where 
$\epsilon _{\theta }$ is the scalar fluctuation dissipation rate. From dimensional analysis, we find that 
$T_{\theta } \propto \epsilon _{\theta } k^{-1}$. For a statistically stationary state, we can rearrange the terms in (3.1) to obtain a 
$E_{\theta } \propto k^{-3}$ scaling. This derivation is similar to Lance & Bataille (Reference Lance and Bataille1991) though they neglected the net local transfer term and performed dimensional analysis on the velocity fluctuation production by the bubbles instead. Finally, as noted above, the original 
$-5/3$ plateau found in figure 2(b) becomes only a peak for increasing 
$\alpha$, i.e. a pure 
$-5/3$ scaling disappears. For high 
$\alpha$, even when 
$-3$ scaling has developed, the mixing of the bubbles may lead to an overall more homogeneous temperature field which leads to a smaller temperature fluctuation at the scales also before the 
$-3$ scaling, leading to a steeper slope also before that. As we also show in the next section, the onset frequency of 
$-3$ scaling becomes lower for higher 
$\alpha$. If 
$\alpha$ further increases, this leads us to expect that such steeper slope (steeper than the originally 
$-5/3$) before the developed 
$-3$ scaling will be further steepened, closer to 
$-3$.
4. Transition frequencies from 
$-5/3$ to $-3$ scaling
The frequencies corresponding to the onset of the 
$-3$ scaling of the scalar spectra are now examined. When 
$\alpha$ increases from 0 % to 5.2 %, the 
$-5/3$ scaling becomes less pronounced and, when 
$\alpha$ is large enough, it is followed by a 
$-3$ scaling, see figure 3(a) for examples. We identify the onset frequency 
$f_t$ by using the parameterisation
\begin{equation} P_{\theta}(f)=\frac{(f/f_L)^{\zeta_{b}}}{[1+ (f/f_t)^2]^{({\zeta_{b}-\zeta_{a}})/{2}}}, \end{equation}
which captures a transition from one scaling to another. Here 
$f_L$ is a fitting parameter which reflects the height of the curve and 
$f_t$ is the transition frequency from the scaling 
$f^{\zeta _b}$ (for 
$f \ll f_t$) to 
$f^{\zeta _a}$ (for 
$f \gg f_t$) with scaling exponents 
$\zeta _b$ and 
$\zeta _a$, respectively, see figure 3. Here we do not fit these exponents, but set 
$\zeta _b = -5/3$ (Kolmogorov–Obukhov value) (Monin & Yaglom Reference Monin and Yaglom1975) and 
$\zeta _a = -3$, consistent with the exponent of the velocity spectra in high-
$Re$ bubbly flows (Risso Reference Risso2018) and reflecting the observed limiting cases. To determine the range that covers most of the 
$-5/3$ and 
$-3$ scaling, we examine the local logarithmic slope 
$\xi (f)$ of the scalar power spectra, which is given by
\begin{equation} \xi(f) = \frac{{\rm d} \log_{10} P_{\theta}(f) }{{\rm d} \log_{10} f}. \end{equation} To estimate the local logarithmic slope at different frequency, for each spectrum, a moving fit of a straight line in logarithmic space (centred at a frequency with a window size of one-quarter decade of the frequency) is employed in order to obtain a local logarithmic slope at that frequency. The results of 
$\xi (f)$ for different 
$\alpha$ are shown in figure 3(b). In the figure, we see that around 1.2 Hz, 
$-5/3$ scaling starts to develop and the 
$-3$ scaling ends around 15 Hz for 
$\alpha = 2.2\,\%$. Therefore, 1.2–15 Hz is a reasonable range that covers the transition from 
$-5/3$ to 
$-3$ scaling for all cases.
We fit the spectra of figure 2(b) with (4.1) with 
$\zeta _b=-5/3$ and 
$\zeta _a = -3$ only for the 
$\alpha \geq 1.0\,\%$ cases, because for 
$\alpha < 1\,\%$ the transition frequencies resulting from the fits are either outside the fitting range or very close to the boundary of the fitting range. The fits for 
$\alpha = 1.0\,\%$ and 
$\alpha = 5.2\,\%$ are shown in figure 3(a). It shows qualitatively nice fits for 
$\alpha = 1.0\,\%$, even though this case does not exhibit a fully developed 
$-3$ scaling. The figure also shows that for the case of 
$\alpha = 5.2\,\%$ the scaling is close to 
$-3$ after the transition and that the transition frequency 
$f_t$ is smaller for this higher value of 
$\alpha$.
Table 2 lists the transition frequencies 
$f_t$ obtained from the parameterisation (4.1) as a function of 
$\alpha$ for 
$\alpha \geq 1.0\,\%$. A decreasing trend of 
$f_t$ with increasing 
$\alpha$ is found. Currently, there is no prediction for the onset frequencies 
$f_t$ of the bubble-induced subrange for the passive scalar spectra. However, there is a candidate of the onset frequency of the 
$-3$ scaling in the energy spectra, which is proportional to 
$V_r/d_{eq,V}$, where 
$V_r$ is the relative (to the liquid phase) rise velocity of the bubbles and 
$d_{eq,V}$ is the mean volume equivalent diameter of the bubbles (Alméras et al. Reference Alméras, Mathai, Lohse and Sun2017). This frequency is given by 
$f_{{bub}} \equiv V_r/\lambda _{{bub}} \equiv C_D V_r/d_{eq,V}$ (Alméras et al. Reference Alméras, Mathai, Lohse and Sun2017), where 
$C_D$ is the drag coefficient of the bubbles and 
$\lambda _{{bub}} \equiv d_{eq,V}/ C_D$ is the candidate for the onset scale for the 
$-3$ scaling in the energy spectra in the wavenumber space (Riboux, Risso & Legendre Reference Riboux, Risso and Legendre2010; Risso Reference Risso2018). Apart from 
$f_{{bub}}$, there are several other frequency scales that are proportional to 
$V_r/d_{eq,V}$, which roughly locate a ‘bump’ in the energy spectra just before the 
$-3$ scaling for high enough 
$\alpha$, summarised in Alméras et al. (Reference Alméras, Mathai, Lohse and Sun2017). For example, based on their observations in their energy spectra, Alméras et al. (Reference Alméras, Mathai, Lohse and Sun2017) gives
\begin{equation} f_{c,A} = 0.14 \langle V_r \rangle /\langle d \rangle, \end{equation}
where we stress again that 
$d$ is the area-equivalent diameter whereas in Riboux, Legendre & Risso (Reference Riboux, Legendre and Risso2013) the mean volume-equivalent diameter is used instead. Table 2 includes these values of 
$f_{c,A}$, which we are also shown in figure 3(c). As the 
$-3$ scaling has not developed for low 
$\alpha$ and, thus, the identification of 
$f_t$ for low 
$\alpha$ is hard to determine, we can only conclude that both 
$f_c$ and 
$f_t$ identified here have decreasing trends but of the same order of magnitude in our parameter regime.
Apart from 
$f_{c,A}$, when 
$\alpha$ becomes large, the typical time scale for the effective diffusion for passive scalar in bubble-induced turbulence is the time scale for successive bubble passages (Alméras et al. Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015, Reference Alméras, Mathai, Sun and Lohse2019), which is given by
\begin{equation} T_{2B} = \frac{d}{\alpha V_r \chi^{2/3}}. \end{equation}
The numeral values of the quantities on the right-hand side of the above equation can be found in table 2. For 
$\alpha =4.7\,\%$, 
$T_{2B} = 0.39$ s which results in a frequency of 
$2.6 \pm 0.7$ Hz, which is half the size than 
$f_{c,A} = 4.9 \pm 0.5$ Hz. The transition frequency 
$f_t = 3.9 \pm 0.2$ Hz for the corresponding 
$\alpha$. For low void fraction (
$\alpha =0.6\,\%$), we find 
$T_{2B} = 1.57$ s, giving a frequency of 
$0.7 \pm 0.3$ Hz. Therefore, for high 
$\alpha$, this is comparable to our 
$f_t$. However, we see that the trend of this prediction is going up with increasing 
$\alpha$, whereas we find that the transition frequency coming from the spectra goes down with 
$\alpha$ (see figure 3c). This is consistent with that for large 
$\alpha$ the time scale for effective diffusion by bubble-induced turbulence can be dominated by the bubble passage time; whereas for small 
$\alpha$ such effective diffusive time scale can be, for example, the Lagrangian integral time scale 
$T_L$ (Alméras et al. Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015, Reference Alméras, Mathai, Sun and Lohse2019). The cross-over of the relevant time scale can occur when 
$T_{{2B}}= T_L$ (Alméras et al. Reference Alméras, Mathai, Sun and Lohse2019). However, the Lagrangian integral time scales for 
$\alpha >0$ is not available in our study. We cannot exclude the time scales mentioned previously as relevant to the onset of the 
$-3$ scaling in the thermal spectra in bubbly flows.
5. Conclusion
The speculation on the 
$-3$ scaling mechanism for a scalar spectrum requires diffusivity to be important. This means that the inverse of scalar diffusive time scale needs to be around the onset frequency 
$f_t$ of the 
$-3$ subrange. A direct verification of this hypothesis may be accomplished by investigation of the diffusive and spectral transfer budgets of the scalar fluctuations at the 
$-3$ subrange by direct numerical simulations similar to what (Pandey et al. Reference Pandey, Ramadugu and Perlekar2020) did. Our hypothesis will be confirmed if within the 
$-$3 subrange, the spectral transfer of the scalar fluctuation (behaving as 
$T_{\theta } \propto k^{-1}$) is balanced by the diffusive dissipation (
$2\kappa k^2E_{\theta }(k)$).
To conclude, when a swarm of high-Reynolds-number bubbles are injected into a turbulent thermal mixing layer (
$0.3\,\% \leq \alpha \leq 5.2\,\%$), the scalar spectra that originally have a 
$-5/3$ scaling, for large frequencies develop a steeper 
$-3$ slope for larger gas concentration 
$\alpha$. As opposed to the energy spectra, which abruptly developed a 
$-3$ scaling with the enhancement of the small-scale content for non-zero 
$\alpha$, the change of spectral scaling for passive scalar spectra is more gradual by smoothing out the small-scale content when 
$\alpha$ increases. Although the 
$-3$ scaling at 
${O}(100 \text {Hz})$ of the energy spectra is attributed to the balance of the production of velocity fluctuations due to bubbles and the molecular viscous dissipation, we speculate that the physical mechanism of the steeper spectral slope in the passive scalar spectra is due to enhanced mixing of bubbles that promotes the homogenisation of small-scale temperature differences. The scaling saturates to 
$-3$ around 
${O}$(10 Hz) when the smoothing of temperature fluctuations is so strong such that the local net transfer of the spectral fluctuation is directly diffused by molecular diffusivity. The transition frequency 
$f_t$ from the 
$-5/3$ to 
$-3$ scaling in the scalar spectra are obtained from the parameterisation (4.1) and found to monotonically decrease with increasing 
$\alpha$.
In our experiments we had 
${Pr} = 6.5$ at 
$22.7\,^\circ$C. This 
$Pr$ could experimentally be changed in a certain narrow regime by changing the water temperature (e.g. 
${Pr} = 12.9$ and 
${Pr} = 2.2$). To explore the 
$Pr$ dependence, direct numerical simulations would be most welcome. Note however that the much smaller 
$Pr$ numbers of liquid metal (e.g. 
${Pr} \approx 0.006$ for liquid sodium at 
$570$ K; see Cioni et al. (Reference Cioni, Horanyi, Krebs and Müller1997)) are again experimentally accessible and that the fluctuations of the temperature in turbulent bubbly sodium may be relevant for applications in working facilities in nuclear power stations.
Acknowledgements
We thank D.P.M. van Gils, G.-W. Bruggert and M. Bos for technical support. We thank T. Chan and J. Cheung for preliminary measurements.
Funding
This work was supported by The Netherlands Center for Multiscale Catalytic Energy Conversion (MCEC), an NWO Gravitation Programme funded by the Ministry of Education, Culture and Science of the Government of The Netherlands. Chao Sun acknowledges the financial support from Natural Science Foundation of China under grant nos. 11988102 and 91852202.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Method of calculating the energy and scalar spectra in two-phase flow
The hot-film voltage 
$E_{hf}(t)$ is measured over time. We remove parts of the signal where bubbles interact with the sensor using a threshold method on the time derivative of the measured voltage across the hot-film (see Rensen, Luther & Lohse Reference Rensen, Luther and Lohse2005a and the references therein) or the time derivative of the liquid velocity (Alméras et al. Reference Alméras, Mathai, Lohse and Sun2017). During a bubble collision with a probe, there are rapid changes of 
$E_{hf}$, and the typical shapes of 
$E_{hf}(t)$ when bubbles colliding with a hot film can be found in, for example, Rensen et al. (Reference Rensen, Luther and Lohse2005a,Reference Rensen, Luther, de Vries and Lohseb). In this work, by inspection of the derivative 
${\rm d}E_{hf}/{\rm d}t$ over time, a fixed threshold 
$|{\rm d}E_{hf}/{\rm d}t| \geq 500$ V s
$^{-1}$ is used to detect the moments when bubbles are colliding the probe. We find that the findings and the conclusions of this work are robust against the threshold values within a reasonable range. Apart from the said threshold, we also set the minimal time between two successive bubble collisions (5 ms) in order to determine the moments when bubbles impinging and leaving the probe and to capture the entire bubble interaction events. We check the latter threshold by examining the distributions of the bubble residence time and any unphysical bumps present in the energy spectra due to many unphysically short bubble collisions. To calculate the velocity power spectra, we use the Bartlett method (Oppenheim & Schafer Reference Oppenheim and Schafer2009) with linear interpolations between the ‘gaps’ of the liquid phase velocity over time for each 
$\alpha$. Each partition segment (20 s) cover at least one turnover time by inspecting the auto-correlation function of the signal. Linear interpolation between the gaps of the liquid signal to calculate the power spectra was also used in Alméras et al. (Reference Alméras, Cazin, Roig, Risso, Augier and Plais2016, Reference Alméras, Mathai, Lohse and Sun2017). The discussion on the effects of using linear interpolation can be found in Rensen et al. (Reference Rensen, Luther and Lohse2005a) and Alméras et al. (Reference Alméras, Cazin, Roig, Risso, Augier and Plais2016). For the optical method that calculated the energy spectra while also having evenly sampled time of the liquid velocity field, recently Ma et al. (Reference Ma, Hessenkemper, Lucas and Bragg2022) employed high-resolution particle shadow velocimetry. The advantage of this method is that the two-dimensional velocity field is captured simultaneously without employing Taylor hypothesis (Ma et al. Reference Ma, Hessenkemper, Lucas and Bragg2022). However, the method in this work only measured single-point velocity field and requires Taylor hypothesis when the time coordinate is converted into spatial coordinate (and it is only one-dimensional).
For the computation of the scalar spectra, we did not create gaps of the temperature time series to distinguish between the gas and liquid phases. This is mainly due to that the response time of the thermistor is comparable to the bubble residence time and the bubbles have similar temperature as the surrounding liquid, which is discussed and shown in chapter 5 in Dung (Reference Dung2021). The Bartlett method (Oppenheim & Schafer Reference Oppenheim and Schafer2009) is applied to the temperature signal using partition segments of 10 s which cover at least one turnover time which was found by inspecting the auto-correlation function of the signal.
Appendix B. Obtaining the thermal integral length scale
For the thermal integral length scale, intuitively it is of the order of the size of the thermal mixing layer. For a single-phase thermal-mixing layer, the mean temperature profile follows a self-similar shape that can be characterised by a complementary error function (Ma & Warhaft Reference Ma and Warhaft1986; de Bruyn Kops & Riley Reference de Bruyn Kops and Riley2000). Therefore, quantitatively we first fit the mean temperature profile (with an offset of the mean temperature at the coldest position) 
$\tilde {T} (x)$ by the function
\begin{equation} \tilde{T} (x)= \frac{{\rm \Delta} T}{2}\text{erfc} \left(\frac{x-x_c}{L_{ \theta}/2}\right), \end{equation}
where 
$x_c$ is the 
$x$-coordinate at the centre of the measurement section, 
$\text {erfc}(x)$ is the complementary error function; and 
${\rm \Delta} T$ and 
$L_{\theta }/2$ are the fitting parameters which describe the asymptotic value of the mean profile and the characteristic length of the thermal mixing layer, respectively, see figure 4 for the fit. Then 
$L_{\theta }$ is defined as the integral length scale of the scalar fluctuations. The results are summarised in table 1.
The mean temperature profile 
$\tilde {T}$ versus 
$({x-x_c})/{L_{\theta }/2}$ for 
$\alpha = 0\,\%$ and the corresponding fit using (B1) (solid black line), where 
$\tilde {T} \equiv \langle T - T_{{inlet}} \rangle -\langle T - T_{{inlet}} \rangle |_{x/L = 0.9}$, 
$T_{{inlet}}$ is the temperature measured of the incoming liquid before being heated up by the heaters, 
$x_c = L/2$ is the 
$x$-coordinate of the middle of the tunnel and 
$L_{\theta }$ is the fitting parameter that measures the width of the thermal mixing layer (B1).
Open access
