1. Introduction
The transport of scalar fields, such as heat or chemical species, in bubbly flows occurs in many industrial and natural processes such as chemical reactors, heat exchangers and atmosphere–ocean exchanges. Empirical evidence suggests that bubbly flows enhance the transport of scalars without moving parts, diminishing costs (Mudde Reference Mudde2005).
Bubbles rising in a liquid at moderate volume fractions induce peculiar velocity fluctuations commonly known as bubble-induced turbulence (Lance & Bataille Reference Lance and Bataille1991; Risso Reference Risso2018). Experiments (Mercado et al. Reference Mercado, Gomez, Van Gils, Sun and Lohse2010; Riboux, Risso & Legendre Reference Riboux, Risso and Legendre2010; Mendez-Diaz et al. Reference Mendez-Diaz, Serrano-Garcia, Zenit and Hernandez-Cordero2013) and numerical studies (Pandey, Ramadugu & Perlekar Reference Pandey, Ramadugu and Perlekar2020; Innocenti et al. Reference Innocenti, Jaccod, Popinet and Chibbaro2021) have shown a robust power-law scaling with an exponent of 
$-3$ in the kinetic energy spectrum at an approximate interval of wavenumbers 
$k$ below the bubble diameter before viscous dissipation occurs. This scaling differs from the classical 
$-5/3$ observed in single-phase homogeneous isotropic turbulence.
Although the properties of bubble-induced turbulence have been thoroughly investigated, the dynamics and statistics of a passive scalar in such flows have only recently received attention. For example, Alméras et al. (Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015) examined experimentally the mixing of a low-diffusive dye in homogeneous bubbly flows and showed that the scalar dispersion could be modelled as an anisotropic diffusion process with the effective diffusivity 
$\propto \phi ^{0.4}$ at low gas volume fractions 
$\phi$. A follow-up work (Alméras et al. Reference Alméras, Cazin, Roig, Risso, Augier and Plais2016) shows, in the same experimental configuration, a 
$-3$ scaling of the scalar spectrum in the frequency domain. Loisy, Naso & Spelt (Reference Loisy, Naso and Spelt2018) studied numerically the passive scalar mixing in bubbly flows of up to 
$12$ bubbles at low bubble Reynolds numbers (
$Re=30$). That study showed that the convective contribution (due to bubble-induced agitation) to the effective scalar diffusivity is dominant for most common bubbly flows. Gvozdić et al. (Reference Gvozdić, Alméras, Mathai, Zhu, van Gils, Verzicco, Huisman, Sun and Lohse2018) experimentally investigated the heat transport in a bubble column heated on one lateral side and cooled on the other. They found the effective thermal diffusivity in the horizontal direction 
$\propto \phi ^{0.45}$ for 
$\phi \leq 5\,\%$. They also observed a scaling of 
$-1.4$ in the temperature spectrum at frequencies around 
$0.1\unicode{x2013}3$ Hz; however, the thermistors used could not resolve the higher frequencies present in the bubble-induced turbulence. Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023) studied experimentally the thermal spectrum scaling of a thermal mixing layer in vertical channel bubbly flow with an active turbulent grid. They showed a spectrum scaling transition from 
$-5/3$ to a 
$-3$ scaling in the frequency domain for large enough 
$\phi$ that clearly shows the bubbles influence the thermal spectra. The work provides hypotheses and scaling arguments for the existence of the 
$-3$ scaling in the scalar spectrum.
For the first time, we investigate, using multiphase direct numerical simulations (DNS), the spectral scaling of a passive scalar 
$c$ in bubble-induced turbulence. Specifically, we study the influence of the liquid Schmidt number and 
$\phi$ on the scalar spectra (at length scales of 
$O(10)$ bubble diameters down to the viscous dissipation scales) and the effective scalar diffusivity. Statistically steady scalar fluctuations are generated by imposing constant scalar gradients in both the vertical and horizontal directions. Furthermore, we compute the energy budget of the scalar spectra to assess the hypothesis proposed in Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023) that the spectral scalar transfer 
$T_c(k) \propto k^{-1}$ where the scalar spectra show a 
$-3$ scaling. In addition, we compare the results from our bubbly flow simulations with those obtained in single-phase isotropic turbulence to elucidate the effects of the bubbles on the scalar dynamics.
2. Methodology
2.1. Problem statement
We numerically study the statistics of a passive scalar field with an imposed constant gradient in a fully periodic homogeneous bubbly flow domain where the bubbles agitate an initially quiescent liquid. We assume the scalar field to be continuous across the bubble interfaces and the surface tension, densities, viscosities and scalar diffusivities constant in the two phases. All variables are made non-dimensional using the spherical equivalent bubble diameter 
$d_b$, characteristic rise velocity 
$\sqrt {g d_b}$ (
$g$ is the gravitational acceleration) and the liquid density 
$\rho _l$. The problem is completely described by the seven, a priori known, dimensionless parameters; the Eötvös number 
$Eo = {\rho _l g d_b^2}/{\sigma }$ relating buoyancy to surface tension forces, the Galilei number 
$Ga = {\rho _l \sqrt {gd_b}d_b}/{\mu _l}$ that is the ratio of buoyancy to viscous forces, the density ratio 
$\rho _r={\rho _g}/{\rho _l}$, the dynamic viscosity ratio 
$\mu _r = {\mu _g}/{\mu _l}$, the gas volume fraction 
$\phi = N_b{\rm \pi} d_b^3/(6L^3)$ and the Schmidt numbers 
$Sc_l = \nu _l/D_l$ and 
$Sc_g = \nu _g/D_g$. Here, 
$\sigma$ is the surface tension, 
$N_b$ the number of bubbles, 
$L$ the side length of the cubic domain, 
$\nu =\mu /\rho$ the kinematic viscosity and the subscripts 
$l$ and 
$g$ denote the liquid and gas phases. We use 
$N_b \geq 40$ monodisperse bubbles to get statistics independent of 
$N_b$ (Loisy, Naso & Spelt Reference Loisy, Naso and Spelt2017). We choose 
$\rho _r=10^{-3}$ and 
$\mu _r=10^{-2}$ that resemble air–water systems. These ratios are very small and the exact values are physically insignificant for most gas–liquid systems of practical interest (Bunner & Tryggvason Reference Bunner and Tryggvason2002). We fix 
$Ga=390$ and 
$Eo=0.85$ that correspond to 
$2.5$ mm air bubbles in water. These parameters are characteristic of practically relevant systems and several experimental studies (Riboux et al. Reference Riboux, Risso and Legendre2010; Mendez-Diaz et al. Reference Mendez-Diaz, Serrano-Garcia, Zenit and Hernandez-Cordero2013; Alméras et al. Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015; Gvozdić et al. Reference Gvozdić, Alméras, Mathai, Zhu, van Gils, Verzicco, Huisman, Sun and Lohse2018; Dung et al. Reference Dung, Waasdorp, Sun, Lohse and Huisman2023).
We study a scalar field, such as the concentration of a chemical species or the temperature, that can be assumed to be a passive scalar if the effect of the scalar on the fluid properties (such as viscosity and density) is small (and the effects of viscous heating can be ignored) (Gvozdić et al. Reference Gvozdić, Alméras, Mathai, Zhu, van Gils, Verzicco, Huisman, Sun and Lohse2018; Loisy et al. Reference Loisy, Naso and Spelt2018). Although we choose to use typical parameters in the context of mass transfer, the results are relevant for heat transport, given that the assumptions mentioned above hold. We analyse the effects of the scalar gradient direction, liquid Schmidt numbers 
$Sc_l=(0.7, 1.5, 3, 7)$ and the gas volume fractions 
$\phi =1.7\,\%$ and 
$5.2\,\%$ on the scalar spectra and transport properties. The complete set of DNS parameters is shown in table 1. Note that for each case in table 1, we simulate two independent scalar fields with an imposed constant gradient in the vertical 
$\nabla ^v \langle c \rangle$ and horizontal 
$\nabla ^h \langle c \rangle$ direction, respectively.
Table 1. Non-dimensional parameters of the DNS: 
${Re}_d=V_d d_b / \nu _l$ is defined with the drift velocity 
$V_d$, the average relative velocity between the gas and liquid phases; 
$u_{l,{std}}$ and 
$c_{l,{std}}$ are the total standard deviation of the liquid velocity and liquid scalar fluctuations, respectively; 
$\epsilon _{l,{c}}=D_l\langle \boldsymbol {\nabla } c'_l \boldsymbol{\cdot } \boldsymbol {\nabla } c'_l \rangle$ are the scalar dissipation rates in the liquid and 
${Sh} = {\boldsymbol {D}}_{{conv}}/D_{{mol,s}}$ is the Sherwood number. The superscripts 
$v$ and 
$h$ represent values for the cases with an imposed scalar gradient in the vertical 
$\nabla ^v \langle c \rangle$ and horizontal 
$\nabla ^h \langle c \rangle$ directions, respectively. For the single-phase cases S1–4, despite the isotropy, the scalar statistics are shown in the columns for the horizontal direction.
Additionally, we study the scalar dynamics in single-phase isotropic turbulence to compare with that in bubbly flows. We use a fully periodic domain and the same 
$Sc_l$ values as in the bubbly flow cases. Turbulence is maintained by an isotropic stochastic forcing localised at a wavenumber where the bubble-induced velocity spectrum shows a maximum and the dissipation is kept equal to the one of the bubble suspension 
$\epsilon =\phi g V_d$.
2.2. Numerical method
The scalar field is decomposed into 
$c = \boldsymbol {\nabla }\langle c \rangle \boldsymbol{\cdot } \boldsymbol {x} + c'$, where 
$\langle c \rangle$ represents the mean scalar field that we specify as a constant slope linear field (
$\boldsymbol {\nabla } \langle c \rangle =1$) and 
$c'$ is the scalar disturbance due to the bubbles’ motion that we solve numerically. The bubbly suspension is solved using the volume of fluid approach, where the two phases are tracked using a volume fraction field 
$f$ that is equal to 1 in the liquid and 0 in the gas. The governing equations read
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$
$$\begin{gather}\rho\dfrac{D \boldsymbol{u}}{D t} = (\rho -\langle \rho \rangle) \boldsymbol{g} - \boldsymbol{\nabla} p + \dfrac{1}{Ga}\boldsymbol{\nabla} \boldsymbol{\cdot}(2\mu{\boldsymbol{\mathsf{S}}}) + \dfrac{\kappa\delta_S \hat{\boldsymbol{n}}}{Eo} ,\end{gather}$$
$$\begin{gather}\dfrac{\partial f}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u}f) = 0 , \end{gather}$$
$$\begin{gather}\dfrac{\partial c'}{\partial t} + \underbrace{ \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} c'}_{h_c} = \underbrace{\boldsymbol{\nabla} \boldsymbol{\cdot} (D \boldsymbol{\nabla} (\langle c \rangle + c'))}_{d_c} - \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \langle c \rangle , \end{gather}$$
where 
$\boldsymbol {u}$ is the velocity, 
$\boldsymbol {g}$ the gravity vector, 
$p$ the pressure, 
${\boldsymbol{\mathsf{S}}}=(\boldsymbol {\nabla }\boldsymbol {u} + \boldsymbol {\nabla }\boldsymbol {u}^T)/2$ the rate of deformation tensor and 
$\kappa$ and 
$\hat {\boldsymbol {n}}$ are the interface curvature and normal. The additional body force 
$\langle \rho \rangle \boldsymbol {g}$ in (2.2) prevents the flow from accelerating in the gravitational direction (Bunner & Tryggvason Reference Bunner and Tryggvason2002). The density 
$\rho$ and diffusivity 
$D$ are the arithmetic means with 
$f$ of their single-phase counterparts, while 
$\mu$ is the harmonic mean that better approximates gas–liquid interfaces with continuous shear stress (Tryggvason, Scardovelli & Zaleski Reference Tryggvason, Scardovelli and Zaleski2011).
The governing equations are solved with the open-source code Basilisk (Popinet Reference Popinet2015) in a cubic periodic domain with the side length 
$L/d_b = 10.72$. We discretise the domain with equidistant grid points in each direction. The number of grid points is selected to resolve the smallest length scales of the velocity and scalar fields. The Kolmogorov length scale 
$\eta _u = (\nu ^3/\epsilon )^{1/4}$ is estimated using that the dissipation rate equals the power of the buoyancy force per unit mass 
$\epsilon = \phi g V_d$ (Risso Reference Risso2018). For cases B, with 
$\phi =5.2\,\%$, we obtain 
$\eta _u / d_b=0.021$. Using the DNS resolution criterion reported in Pope (Reference Pope2001), a grid spacing 
$\varDelta$ satisfying 
$\varDelta /\eta _u \lessapprox 2.1$ gives a good resolution of the smallest turbulent scales. These values indicate that a grid spacing of about 
$24 \varDelta / d_b$ is sufficient for resolving the velocity field.
The Batchelor length scale is estimated according to 
$\eta _c = \eta _u / Sc_l^{1/2}$ (Batchelor Reference Batchelor1959). This relation indicates that, for the cases B with 
$Sc_l \approx 1$, the resolution of 
$24 \varDelta / d_b$ is sufficient, while for the cases with 
$Sc_l = 3$ and 
$7$, we need approximately 
$40$ and 
$60 \varDelta / d_b$, respectively. The same analysis for cases A, with 
$\phi =1.7\,\%$, shows that 
$Sc_l = 3$ and 
$7$ require about 
$30$ and 
$50 \varDelta / d_b$, respectively.
However, the DNS resolution criterion by Pope (Reference Pope2001) is based on studies on single-phase isotropic turbulence, and, in principle, even smaller scales may be present in the bubbly flow velocity and scalar fields (such as thin boundary layers). For these reasons (and because we evolve multiple scalar fields in a single DNS), we ran all cases with 
$Sc_l = 0.7$ and 
$1.5$ with 
$48 \varDelta / d_b$ and the cases with 
$Sc_l = 3$ and 
$7$ using another refinement level corresponding to 
$96 \varDelta / d_b$.
The Basilisk code has been validated and used extensively for multiphase DNS of bubbly flows (Innocenti et al. Reference Innocenti, Jaccod, Popinet and Chibbaro2021; Hidman et al. Reference Hidman, Ström, Sasic and Sardina2022). The code features a finite volume solver using a time-splitting projection method with standard second-order gradient discretisation and the Bell–Colella–Glaz second-order upwind scheme for the velocity and scalar advection. We use the cell-centred Cartesian multigrid solver that efficiently solves the Helmholtz–Poisson-type problems for the velocity components and the Poisson equation for the pressure correction (Popinet Reference Popinet2009). This allows us to perform high-resolution three-dimensional simulations with a uniform grid at a feasible computational cost. The surface tension force in (2.2) is discretised with a well-balanced method, and the height-function approach is used to compute the interface curvature (Popinet Reference Popinet2018). The piecewise-linear interface reconstruction technique is used to advect the volume fraction field without interface smearing (Scardovelli & Zaleski Reference Scardovelli and Zaleski1999). To avoid bubble coalescence, we implement a repulsive force by locally increasing the surface tension to 
$\sigma _{{rep}} = 2.1\sigma$ (Talley, Zimmer & Bolotnov Reference Talley, Zimmer and Bolotnov2017) only at the part of the bubble interface that is less than 
$d_b$ from another bubble's centre of mass.
The single-phase simulations are performed with a pseudo-spectral solver (Sardina et al. Reference Sardina, Picano, Brandt and Caballero2015). The grid increases from 384 to 640 collocation points with increasing 
$Sc$.
3. Results
This section presents the passive scalar spectra and statistics from our DNS of single-phase isotropic turbulence and homogeneous bubbly turbulent suspensions. The simulation cases and parameters are shown in table 1.
3.1. Characteristics of the scalar dynamics
Figure 1 shows instantaneous snapshots from case B4 of the full periodic computational domain and the bubble interfaces coloured by the total scalar field 
$c$ with 
$\nabla ^v \langle c \rangle$ in figure 1(a) and 
$\nabla ^h \langle c \rangle$ in figure 1(b). Vertical cross-sectional views of the total scalar field 
$c$ and the scalar disturbance field 
$c'$ from figure 1(a) are shown in figure 2 and from figure 1(b) in figure 3.
Figure 1. Istantaneous snapshot from case B4 with bubble interfaces and periodic domain boundaries coloured by the total scalar field 
$c$. The volume fraction is 
$\phi =5.2\,\%, Sc_l = 7$ and 
$\nabla ^v \langle c \rangle =1$ is imposed in panel (a) and 
$\nabla ^h \langle c \rangle =1$ in panel (b).
Figure 2. Vertical cross-sectional view of the total 
$c$ and disturbance 
$c'$ scalar fields from case B4 with 
$\phi =5.2\,\%, \nabla ^v \langle c \rangle =1$ and 
$Sc_l = 7$. The thin black lines represent the bubble interfaces.
Figure 3. Vertical cross-sectional view of the total 
$c$ and disturbance 
$c'$ scalar fields from case B4 with 
$\phi =5.2\,\%, \nabla ^h \langle c \rangle =1$ and 
$Sc_l = 7$.
The generation of the scalar disturbances stems from the last term on the right-hand side of (2.4) that (given a positive 
$\nabla ^v \langle c \rangle$) predicts large negative disturbances in high upward velocity regions and positive disturbances in downward velocity regions. Consequently, we speculate that the scalar disturbances in the case of 
$\nabla ^v \langle c \rangle$ (shown in figure 2b) are generated by several interacting mechanisms. The first is the generation by the gas phase that preferentially moves in the vertical direction (generation at scales 
$k/k_{d_b}\approx 1$). The second is the associated generation by the upward flowing liquid in the bubble wakes that is related to the capture–release mechanism described in Alméras et al. (Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015, Reference Alméras, Risso, Roig, Plais and Augier2018) (generation at scales 
$k/k_{d_b} \gtrapprox 1$). The third is the generation by the average downward flowing liquid in between the rising bubbles due to continuity (generation at scales 
$k/k_{d_b} \lessapprox 1$). The last is the production by the bubble-induced turbulence that was found to dominate the dispersion of a low-diffusive dye in a homogeneous bubbly flow similar to the flow considered in the present study (Alméras et al. Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015).
In contrast, the scalar disturbances in the case of 
$\nabla ^h \langle c \rangle$ (shown in figure 3b) show different characteristics. Here, the scalar disturbances are mainly generated by the bubble-induced turbulence since the average bubble motion in the horizontal directions equals zero. However, since the bubble-induced turbulence is characterised by anisotropy (due to the preferential motion of the bubble in the vertical direction), the influence of the bubble-induced turbulence on the scalar dynamics is not, in general, the same in cases with different mean gradient orientations. Given the different characteristics of the scalar disturbances in the cases of 
$\nabla ^v \langle c \rangle$ and 
$\nabla ^h \langle c \rangle$, we thus expect some differences in the scalar statistics that are analysed in § 3.4.
3.2. Bubble clustering
Following the methodology of Tagawa et al. (Reference Tagawa, Roghair, Prakash, van Sint Annaland, Kuipers, Sun and Lohse2013) and Pandey et al. (Reference Pandey, Ramadugu and Perlekar2020), we analyse the bubble clustering in our DNS with volume fractions 
$\phi =1.7\,\%$ and 
$\phi =5.2\,\%$. We compute the centre of mass of each bubble in the DNS at a statistically steady state and compute Voronoi tessellations using the Voro++ library (Rycroft Reference Rycroft2009). For each volume fraction, we compute the standard deviation 
$\varSigma$ of the Voronoi volumes. This 
$\varSigma$ is compared with the standard deviation 
$\varSigma _{{rnd}}$ of the Voronoi volumes in 200 configurations of randomly positioned, non-overlapping bubbles in the same cubic domain with length 
$L$ as used in the DNS. Tagawa et al. (Reference Tagawa, Roghair, Prakash, van Sint Annaland, Kuipers, Sun and Lohse2013) shows that the clustering indicator 
$\mathcal {C}=\varSigma /\varSigma _{{rnd}}$ quantitatively identifies different clustering morphologies. A value of 
$\mathcal {C}<1$ indicates a regular lattice arrangement, 
$\mathcal {C}=1$ a random bubble distribution and 
$\mathcal {C}>1$ irregular clustering. For our case with 
$\phi =1.7\,\%$, we obtain 
$\mathcal {C}=0.95$ and for 
$\phi =5.2\,\%$ we have 
$\mathcal {C}=1.2$. These values indicate random or weakly irregular clustering in our DNS (Pandey et al. Reference Pandey, Ramadugu and Perlekar2020). Consequently, we observe no significant effect of clustering on our results.
3.3. Velocity spectra
We define the velocity spectrum 
$E_u$ of the suspension at wavenumber 
$k=\lvert {\boldsymbol k}\rvert$ as
\begin{equation} E_u(k) =\frac{1}{2} \sum_{k-{\rm \Delta} k/2 < \lvert{\boldsymbol k}\rvert < k+{\rm \Delta} k/2} \dfrac{\langle \hat{\boldsymbol{u}}'({\boldsymbol k}) \boldsymbol{\cdot} \hat{\boldsymbol{u}}'^*{(\boldsymbol k)} \rangle}{{\rm \Delta} k}, \end{equation}
where ‘
$\hat { }$’ represents the Fourier mode, ‘
$^*$’ is the complex conjugate and 
$\langle {\cdot } \rangle$ is the ensemble average. Similarly, we can define the velocity spectrum of the liquid phase considering the liquid velocity 
${\boldsymbol {\tilde {u}_l}=\boldsymbol {u}\lvert {(\,f=1)}}$; however, we regularise the liquid velocity directly with the volume fraction field to avoid discontinuities that can lead to possible Gibbs phenomena in the Fourier transform defining the regularised liquid velocity as 
$\boldsymbol {u_l}=f\boldsymbol {u}$
\begin{equation} {E}_{u,l}(k) =\frac{1}{2} \sum_{k-{\rm \Delta} k/2 < \lvert{\boldsymbol k}\rvert < k+{\rm \Delta} k/2}\dfrac{\langle \widehat{\boldsymbol{u_l}}'({\boldsymbol k}) \boldsymbol{\cdot} \widehat{\boldsymbol{u_l}}'^*{(\boldsymbol k)} \rangle}{{\rm \Delta} k}.\end{equation} Figure 4(a) shows the suspension and liquid velocity spectra for our cases, and previous bubbly flow DNS by Pandey et al. (Reference Pandey, Ramadugu and Perlekar2020) at similar governing parameters (case R7 in that study). All cases show a peak at scales close to the bubble diameter (
$k/k_{d_b}\approx 1$, with 
$k_{d_b}=2{\rm \pi} /d_b$) whereas at 
$k/k_{d_b} > 1$, the velocity spectra of the bubbly flows scale approximately as 
$k^{-3}$, in agreement with previous studies. The single-phase velocity spectra show an even steeper slope since, here, the velocity fluctuations are only produced at the forcing length scale. This is opposed to bubbly flows where velocity fluctuations are produced by both large scales, and at 
$k/k_{d_b} > 1$ where the fluctuations are continuously produced and directly dissipated in the bubble wakes (Lance & Bataille Reference Lance and Bataille1991). The single-phase velocity spectrum does not show the Kolmogorov 
$-$5/3 scaling in a large interval since the Taylor–Reynolds number is only 39. The 
$-3$ scaling in the bubbly flows is more pronounced in the liquid than in the suspension, so we focus on the scalar dynamics in the liquid phase.
Figure 4. Normalised velocity spectra 
$E_u(k)$ and scalar spectra 
$E_{c,l}(k)$ against the normalised wavenumber 
$k/k_{d_b}$. (a) Velocity spectra of single-phase, bubble suspensions and liquid phases at two different volume fractions and comparisons with Pandey et al. (Reference Pandey, Ramadugu and Perlekar2020). (b) Value of 
$E_c(k)$ for the single phase isotropic turbulence case. (c) Bubbly flow cases with 
$\nabla ^v \langle c \rangle$. (d) Bubbly flow cases with 
$\nabla ^h \langle c \rangle$. The gas volume fraction is 
$\phi =1.7\,\%$ in all cases A1–A5. Case B4 with 
$\phi =5.2\,\%$ is included to illustrate the effects of 
$\phi$ on 
$E_{c,l}(k)$.
3.4. Scalar spectra
We define the suspension and liquid scalar spectra similar to (3.1) and (3.2)
\begin{align} E_c(k) = \sum_{k-{\rm \Delta} k/2 < \lvert{\boldsymbol k}\rvert < k+{\rm \Delta} k/2} \dfrac{\langle \hat{{c}}'({\boldsymbol k}) \hat{{c}}'^*{(\boldsymbol k)} \rangle}{{\rm \Delta} k}, \quad E_{c,l}(k) = \sum_{k-{\rm \Delta} k/2 < \lvert{\boldsymbol k}\rvert < k+{\rm \Delta} k/2} \dfrac{\langle \widehat{{c_l}}'({\boldsymbol k}) \widehat{{c_l}}'^*{(\boldsymbol k)} \rangle}{{\rm \Delta} k}. \end{align} The single-phase scalar spectra are shown in figure 4(b). Below the forcing scale, the classical 
$-5/3$ scaling emerges, especially for the highest 
$Sc$. For lower 
$Sc$, the diffusion significantly influences the dynamics below the forcing scale, like in the velocity spectra, decreasing exponentially due to viscous dissipation (Pope Reference Pope2001).
The normalised liquid scalar spectra for the bubbly flow simulations are shown in figure 4(c) for the cases with 
$\nabla ^v \langle c \rangle$ and in figure 4(d) for the cases with 
$\nabla ^h \langle c \rangle$. For brevity, we only show the spectra for all cases at 
$\phi =1.7\,\%$ since the general trends are the same at 
$\phi =5.2\,\%$. Case B4 is, however, included to illustrate the minor effects of a higher 
$\phi =5.2\,\%$ on 
$E_c/(c^2_{{std}}d_b)$. We do not see the different scaling laws at different volume fractions as in the experiments by Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023). However, in that experiments, the liquid turbulence results from the interaction between the bubble-induced turbulence and the incident turbulence generated by an active grid. This experimental system corresponds to a finite value of the bubblance parameter (ratio of the kinetic energy produced by the bubbles to the turbulent kinetic energy in the absence of bubbles, Rensen, Luther & Lohse Reference Rensen, Luther and Lohse2005). This is in contrast with the present DNS results with only bubble-induced turbulence and an infinite bubblance parameter. Indeed, Prakash et al. (Reference Prakash, Mercado, van Wijngaarden, Mancilla, Tagawa, Lohse and Sun2016) showed experimentally how the 
$-3$ scaling of the velocity spectra emerged as the bubblance parameter was modified. This difference between the present DNS and the experiment by Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023) is a possible reason why we do not observe different scaling laws at different volume fractions. In addition, the latter experimental study considers a thermal mixing layer with a nonlinear mean gradient, whereas the present study imposes a constant mean scalar gradient. The influence of a nonlinear gradient on the scalar spectrum is not clear but may also explain the different observations.
In figures 4(c) and 4(d), we observe a maximum of the scalar spectra at the lowest wavenumber, a signature that the scalar is directly forced by the mean gradients. Consistent with the results of Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023), the liquid scalar spectra show two different scalings. The larger scales show a standard 
$-5/3$ scaling only in a narrow range of wavenumbers. However, after a transition scale, the spectrum scaling changes towards the 
$-$3 scaling or lower for all the Schmidt numbers. When the gradient is in the horizontal direction (figure 4d), we observe a faster decay of the scalar spectrum below the bubble diameter wavenumber. The transition wavenumber increases with the Schmidt number, depending on the molecular diffusivity. It coincides approximately with the bubble wavenumber for the smaller 
$Sc$ since, at 
$Sc$ of order one, we expect that the dynamics of the passive scalar transport is similar to the momentum transport. Case A5 is included to assess the effects of increasing 
$Sc_g$ on the scalar dynamics while keeping 
$Sc_l=0.7$ as in case A1; this case corresponds to the same scalar diffusivity in the gas and liquid phases. We observe that also, in this case, the normalised spectrum shows a 
$-$3 scaling after the bubble diameter wavenumber. Figures 4(c) and 4(d) show only minor differences between cases A1 and A5, indicating that 
$Sc_g$ does not significantly influence the liquid scalar fluctuations. This is, however, not true for the statistics of the suspension, as discussed later.
Figure 5 shows the compensated scalar spectra to more clearly illustrate the scalar spectrum scalings and transitions. Figure 5(a) shows that the classical 
$-5/3$ scaling in the single-phase cases emerges in a narrow range of wavenumbers (where the range increase with 
$Sc_l$) after the forcing length scale. This limited range is due to an insufficient Taylor–Péclet number. Defining the Taylor length scale for the scalar field as 
$\lambda _c = \sqrt {6D_l/\varepsilon _c}c_{l,{std}}$ we obtain Taylor–Péclet numbers 
${Pe}_\lambda = u_{l,{std}} \lambda _c / D_l$ in the range of 33–143 and Taylor length scales 
$\lambda _c$ in the range of 0.15–0.6
$d_b$. These values indicate that molecular diffusion significantly affects the scalar dynamics at scales comparable to or slightly smaller than the bubble size 
$d_b$. It is thus reasonable that the inertial range with the 
$-5/3$ scaling is not very pronounced. The 
$-5/3$ scaling is, however, well established for single-phase DNS with wider inertial ranges (Corrsin Reference Corrsin1951; Gotoh & Watanabe Reference Gotoh and Watanabe2012). Since the forcing length scale in the single phase cases is close to the bubble size, it is reasonable that the 
$-5/3$-scaling is not pronounced in the bubbly flow cases either, as shown in the compensated spectra of figure 5(b).
Figure 5. Compensated and normalised scalar spectra 
$E_{c,l}(k)$ against the normalised wavenumber 
$k/k_{d_b}$. (a) Single-phase isotropic turbulence cases compensated by 
$k^{5/3}$. (b) Bubbly flow cases with 
$\nabla ^h \langle c \rangle$ compensated by 
$k^{5/3}$. (c) Bubbly flow cases with 
$\nabla ^v \langle c \rangle$ compensated by 
$k^{3}$.
Figure 5(c) shows the compensated scalar spectra for the bubble flow cases with 
$\nabla ^v \langle c \rangle$. Here, it is clear that the spectra change scaling towards approximately 
$-3$ at scales 
$k/k_{d_b} \geq 1$ and that the transition occurs at higher 
$k$ for increasing 
$Sc_l$.
3.5. Scalar spectral budget
The transition to the 
$-$3 scaling is clearly a footprint of the bubble-induced agitation mechanisms that influence the scalar dynamics also at scales below the bubble diameter 
$k/k_{d_b}\gg 1$. To investigate the origin of this scaling, we compute the spectral budget of the passive scalar by manipulating (2.4) in Fourier space (Monin & Yaglom Reference Monin and Yaglom1975)
\begin{equation} \dfrac{\partial E_c(k,t)}{\partial t} = T_c(k,t) - \mathcal{E}_c(k,t) + P_c(k,t),\end{equation}
where 
$T_c(k,t) = \sum _{(k-{\rm \Delta} k/2 < \lvert {\boldsymbol k}\rvert < k+{\rm \Delta} k/2)} -\langle \widehat {h_c}\hat {c}'^* + \widehat {h_c}^*\hat {c}' \rangle / {\rm \Delta} k$ is the local transfer term, 
$\mathcal {E}_c(k,t) = \sum _{(k-{\rm \Delta} k/2 < \lvert {\boldsymbol k}\rvert < k+{\rm \Delta} k/2)} \langle \widehat {d_c}\hat {c}'^* + \widehat {d_c}^*\hat {c}' \rangle / {\rm \Delta} k$ is the dissipation term and 
$P_c(k,t) = \sum _{(k-{\rm \Delta} k/2 < \lvert {\boldsymbol k}\rvert < k+{\rm \Delta} k/2)} \boldsymbol {\nabla } \langle c \rangle \boldsymbol {\cdot } \langle (\hat {\boldsymbol {u}}'\hat {c}'^* + \hat {\boldsymbol {u}}'^*\hat {c}')\rangle / {\rm \Delta} k$ is the production at wavenumber 
$k$.
Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023) speculated that the passive scalar transport mechanisms are not due to direct scalar production and simultaneous dissipation at scales smaller than the bubble diameter, as for the momentum transport. Instead, scalar production, generated at larger scales by the mean gradient, decays faster so that at small scales just a balance between liquid scalar transfer 
$T_{c,l}$ and spectral dissipation 
$\mathcal {E}_{c,l}(k,t)$ occurs. By dimensional arguments, they show that, if the liquid scalar transfer is only a function of scalar dissipation and wavenumber, then the scaling is a 
$-1$ power law for the wavenumber 
$T_{c,l}=T_{c,l}(\varepsilon _{c,l},k)\propto \varepsilon _{c,l}k^{-1}$. Since the liquid spectral dissipation is defined as 
$\mathcal {E}_{c,l}(k,t)=2D_lk^2E_{c,l}$, it is trivial to show that if a balance exists with the scalar transfer term, then the scalar spectra should scale as a 
$-$3 power law in 
$k, E_{c,l} \propto k^{-3}$. The different terms of (3.4) are challenging to measure in experiments, so, as suggested in Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023), DNS are needed to estimate them. Here, we assess their hypotheses, extracting the relevant statistics from our numerical simulations.
In figure 6(a,b,c), we show the normalised production 
$P_{c,l}(k)/(\epsilon _{c,l} d_b)$ of (3.4) for the isotropic single-phase (a) and bubbly flow cases with 
$\nabla ^v \langle c \rangle$ (b) and 
$\nabla ^h \langle c \rangle$ (c). The panels show that a significant portion of the scalar fluctuations is produced at approximately 
$k/k_{d_b} \leq 1$ for both the single-phase and bubbly flow cases. At 
$k/k_{d_b} > 1, P_{c,l}(k)$ scales as approximately 
$-7/3$ in a narrow range of 
$k$ for the single-phase cases (as predicted by Lumley (Reference Lumley1964) with dimensional analysis and found in the DNS of Gotoh & Watanabe Reference Gotoh and Watanabe2012). Bubble suspensions, instead, show that 
$P_{c,l}(k)$ decays with a faster power law (approximately 
$k^{-3}$). Therefore, at a statistically steady state, (3.4) implies that the scalar spectra 
$E_{c,l}(k)$ at large wavenumbers are governed by the balance of the net transfer and the diffusive dissipation. Interestingly, we observe a self-similarity in the normalised production spectra for all the cases at different liquid diffusivity. This self-similarity can probably be explained by the fact that the production at large scales is due to the mean scalar gradient, while the liquid diffusivity only influences the smaller scales where the production is small.
Figure 6. Normalised scalar production term 
$P_{c,l}(k)/(\epsilon _{c,l} d_b)$ for the single-phase cases (a) and bubbly flow cases A1–A5 with 
$\phi =1.7\,\%$ for 
$\nabla ^v \langle c \rangle$ (b) and 
$\nabla ^h \langle c \rangle$. (c) Normalised scalar transfer term 
$|T_{c,l}(k)|/(\epsilon _{c,l} d_b)$ for the corresponding single-phase cases (d) and bubbly flow cases with 
$\nabla ^v \langle c \rangle$ (e) and 
$\nabla ^h \langle c \rangle$ (f).
The bottom panels of figure 6 show the absolute value of the scalar transfer spectra. When the production becomes significantly smaller than the transfer, the latter scales close to 
$-1$ in a limited range of wavenumbers. In the same range, the scalar spectra show approximately a 
$-3$ scaling. These trends are highlighted in the zoomed-in plots of figure 7 that show the scalar spectrum, absolute transfer and the production spectrum for each case A1–A4 with 
$\nabla ^v \langle c \rangle$. These plots are consistent with the proposed hypothesis in Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023) that, at large wavenumbers (
$k/k_{d_b} \gtrapprox 3$ in figure 7), the scalar spectrum is governed by the balance of the net transfer and the diffusive dissipation (as the production is negligible) and that the 
$-1$ scaling of the transfer term induces the 
$-3$ scaling of the scalar spectra. The trends shown in figure 7 are observed for all our bubble suspension cases.
Figure 7. Zoomed-in plots of the scalar, transfer and production spectra for our cases A1–A4 with 
$\nabla ^v \langle c \rangle$. In the same range where 
$|T_{c,l}(k)| \propto k^{-1}$, the scalar spectra show a transition to approximately 
$E_{c,l} \propto k^{-3}$.
3.6. Effective scalar diffusivity
The effective scalar diffusivity describes the macroscale scalar flux in the suspension seen as a continuum. Following the detailed derivation of Loisy (Reference Loisy2016), the contributions to the effective diffusivity are obtained by ensemble averaging the scalar transport equation
\begin{equation} \dfrac{\partial \langle c \rangle }{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \langle \boldsymbol{q} \rangle = 0 , \end{equation}where the average flux can be formulated as
\begin{align} \langle \boldsymbol{q} \rangle = \langle\boldsymbol{u}\rangle \langle c\rangle - D_{{mol,s}} \boldsymbol{\nabla}\langle c\rangle + (1-\phi)\langle \boldsymbol{u}'c' \rangle_l + \phi\langle \boldsymbol{u}'c' \rangle_g - (1-\phi)D_l \langle \boldsymbol{\nabla} c' \rangle_l -\phi D_g \langle \boldsymbol{\nabla} c' \rangle_g . \end{align}
Here, we have used that the phase-ensemble average of an arbitrary variable 
$A_s$ of the suspension is 
$\langle A_s \rangle = \langle f A_l + (1-f)A_g \rangle = (1-\phi )\langle A_l \rangle + \phi \langle A_g \rangle$. The effective diffusivity tensor is defined as
$$\begin{gather} {\boldsymbol{\mathsf{D}}}_{{eff}} = D_{{mol,s}} {\boldsymbol{\mathsf{I}}} + {\boldsymbol{\mathsf{D}}}_{{conv}} + {\boldsymbol{\mathsf{D}}}_{{diff}}, \end{gather}$$
$$\begin{gather}{\boldsymbol{\mathsf{D}}}_{{conv}} \boldsymbol{\cdot} \boldsymbol{\nabla} \langle c \rangle =-(1-\phi)\langle \boldsymbol{u}'c' \rangle_l - \phi\langle \boldsymbol{u}'c' \rangle_g, \end{gather}$$
$$\begin{gather}{\boldsymbol{\mathsf{D}}}_{{diff}}\boldsymbol{\cdot} \boldsymbol{\nabla} \langle c \rangle = (1-\phi)D_l \langle \boldsymbol{\nabla} c' \rangle_l +\phi D_g \langle \boldsymbol{\nabla} c' \rangle_g, \end{gather}$$so that (3.5) becomes
\begin{equation} \dfrac{\partial \langle c \rangle }{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\langle\boldsymbol{u}\rangle \langle c\rangle) - \boldsymbol{\nabla} \boldsymbol{\cdot} ({\boldsymbol{\mathsf{D}}}_{{eff}} \boldsymbol{\cdot} \boldsymbol{\nabla} \langle c\rangle) =0 . \end{equation} For all cases in this study 
${\boldsymbol{\mathsf{D}}}_{{diff}} < 0.02{\boldsymbol{\mathsf{D}}}_{{conv}}$. The value of 
${\boldsymbol{\mathsf{D}}}_{{conv}}$ is larger also in Loisy et al. (Reference Loisy, Naso and Spelt2018) for lower bubble Reynolds numbers (
$Re\approx 30$) and with the same diffusivity in both phases. We define 
${\boldsymbol {D}}^v_{{conv}}$ as the non-zero diagonal component of 
${\boldsymbol{\mathsf{D}}}_{{conv}}$ for an imposed scalar gradient 
$\nabla ^v \langle c \rangle$ and 
${\boldsymbol {D}}^h_{{conv}}$ for 
$\nabla ^h \langle c \rangle$. All off-diagonal components of 
${\boldsymbol{\mathsf{D}}}_{{conv}}$ are zero (Loisy Reference Loisy2016).
We define the Sherwood number as 
${Sh} = {\boldsymbol {D}}_{{conv}} / D_{{mol,s}}$ that represents the relative importance of the convective contribution to the molecular contribution in the effective diffusivity of the scalar (in heat transfer, the same definition applies to the Nusselt number). In single-phase isotropic turbulence 
${Sh} \propto {Pe}_{{std}} = u_{{std}}l_u/D_{{mol,s}}$ according to theory and DNS (Gotoh & Watanabe Reference Gotoh and Watanabe2012) at high Péclet numbers, where diffusion is negligible (
$l_u$ is the velocity integral length scale). The same scaling is found at high 
${Pe_d=V_d d_b / D_{{mol,s}}}$ (for cases with both 
$\nabla ^v \langle c \rangle$ and 
$\nabla ^h \langle c \rangle$) for bubbly flows without fully developed bubble-induced turbulence and with a constant scalar diffusivity (Loisy et al. Reference Loisy, Naso and Spelt2018). Figure 8(a) shows that, indeed, 
${Sh} \propto {Pe}_{{std}}$ for our single-phase and bubbly flow simulations with 
$\nabla ^h \langle c \rangle$. Here, we have used 
$l_u = d_b$ as the length scale for 
${Pe}_{{std}}$ in the bubbly flow cases and the forcing length scale in the single-phase simulations. The higher 
${Sh}$ for cases with 
$\nabla ^v \langle c \rangle$ (compared with cases with 
$\nabla ^h \langle c \rangle$) can be explained by the preferential vertical motion of the bubbles due to buoyancy. The ratio of the liquid velocity components in the vertical and horizontal directions is 
$u^v_{l,{std}}/u^h_{{l,std}}\approx 1.4$ (where we obtain for 
$\phi =1.7\,\%: u^v_{l,{std}}=0.29, u^h_{l,{std}}=0.21$ and for 
$\phi =5.2\,\%: u^v_{l,{std}}=0.42, u^h_{l,{std}}=0.31$) for all our bubbly flow cases and where the ratios of the vertical to horizontal components are in excellent agreement with experiments (Riboux et al. Reference Riboux, Risso and Legendre2010). Contrarily, in single-phase isotropic turbulence we have 
$u_{l,{std}}=u^v_{l,{std}}=u^h_{l,{std}}=0.3$. These values suggest that the scalar flux in the single-phase turbulence should be similar to or higher than that in the bubbly flow with 
$\nabla ^h \langle c \rangle$ (where the scalar flux is proportional to 
$u^h_{l,{std}}$). The scalar flux in the single-phase turbulence should, however, be less than that in the bubbly flow with 
$\nabla ^v \langle c \rangle$ (where the total scalar flux is governed by 
$u^v_{l,{std}}$ and the additional mechanisms acting in the vertical direction discussed in § 3.1).
Figure 8. Convective contribution to the effective diffusivity for all our simulation cases. The darker blue squares in (a) and (b) indicate the cases A5 and B5 with 
$\nabla ^v \langle c \rangle$ and 
$Sc_g=7$. (a) Value of 
${Sh}$ against 
${Pe}_{{std}}$ based on the total standard deviation of velocity fluctuations in the liquid; (b) 
${Sh}$ against 
${Pe}_{\phi }$ based on a priori known parameters of the suspension. (c) Non-dimensional 
${\boldsymbol {D}}^*_{{conv}} = {\boldsymbol {D}}_{{conv}} / (\sqrt {g d_b}d_b)$ against 
$\phi$ for 
$\nabla ^v \langle c \rangle$ (top panel) and 
$\nabla ^h \langle c \rangle$ (bottom panel). In both panels, the colour of the points represents the same 
${Sc}$ as in figure 6. In the top panel, the proposed scalings for 
${\boldsymbol {D}}^v_{{conv}}$ are the dash-dotted line (lowest 
$D_{{mol,s}}$, case A4) and the dashed line (highest 
$D_{{mol,s}}$, case B1). The dashed line in the bottom panel is the proposed scaling for 
${\boldsymbol {D}}^h_{{conv}}$.
In the bubbly flow cases with 
$\nabla ^v \langle c \rangle$ we observe approximately 
${Sh} \propto {Pe}^{1.15}_{{std}}$, indicating that molecular diffusion is not negligible. Using the scaling 
$u_{{std}} \propto V_0\phi ^{0.4}$ found in Risso & Ellingsen (Reference Risso and Ellingsen2002) and Riboux et al. (Reference Riboux, Risso and Legendre2010) and noting that the rise velocity of a single bubble 
$V_0 / (\sqrt {g d_b})=O(1)$ we can define 
$Pe_{\phi } = \sqrt {g d_b}\phi ^{0.4}d_b/D_{{mol,s}}$ based on a priori known parameters. Figure 8(b) shows the 
${Sh}$ against 
${Pe}_\phi$ for our bubbly flow cases where we again observe 
${Sh}^v \propto {Pe}^{1.15}_{\phi }$ and 
${Sh}^h \propto {Pe}^{1}_{\phi }$. Solving for 
${\boldsymbol {D}}_{{conv}}$ in 
${Sh}$ of the latter scalings gives 
${\boldsymbol {D}}^v_{{conv}} \propto (\sqrt {g d_b}d_b)^{1.15}\phi ^{0.46}/D_{{mol,s}}^{0.15}$ and 
${\boldsymbol {D}}^h_{{conv}} \propto \sqrt {g d_b}\phi ^{0.4}d_b$. These are similar dependencies on 
$\phi$ as found experimentally in Alméras et al. (Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015) for a low-diffusive dye (
${\boldsymbol {D}}_{{conv}}\propto \phi ^{0.4}$ at lower 
$\phi$) and by Gvozdić et al. (Reference Gvozdić, Alméras, Mathai, Zhu, van Gils, Verzicco, Huisman, Sun and Lohse2018) for heat transport (
${\boldsymbol {D}}_{{conv}} \propto \phi ^{0.45}$ up to 
$\phi =5\,\%$).
The value of 
${\boldsymbol {D}}_{{conv}}$ for each simulation case are shown in figure 8(c). In cases A1–A4 and B1–B4 (coloured dots), 
$Sc_l$ is increased from 
$0.7$ to 
$7$ while 
$Sc_g=0.7$ is maintained. For these cases there is a monotonic increase of both 
${\boldsymbol {D}}^v_{{conv}}$ and 
${\boldsymbol {D}}^h_{{conv}}$ with 
$Sc_l$. However, in cases A5 and B5 (grey crosses) we increase 
$Sc_g=7$ while specifying 
$Sc_l=0.7$. These parameters result in the highest 
${\boldsymbol {D}}^v_{{conv}}$ for the respective cases A1–5 and B1–5 but a 
${\boldsymbol {D}}^h_{{conv}}$ closer to lower 
$Sc_l$.
The influence of 
$Sc_g$ on 
${\boldsymbol{\mathsf{D}}}_{{conv}}$ can be explained by considering the bubbles as a source of scalar disturbance 
$c'$. When a bubble moves in the direction of the imposed scalar gradient, the scalar disturbance in the bubble increases proportional to the last term on the right-hand side of (2.4). The time it takes for the scalar in the bubble to reach an equilibrium with the liquid is proportional to the characteristic equalisation time 
$t_e = R^2/D_{g}$. A higher 
$Sc_g$ thus implies a longer 
$t_e$ and therefore a higher average value of 
$c'$ in the gas phase contribution 
$\phi \langle \boldsymbol {u}'c' \rangle _g$ to 
${\boldsymbol{\mathsf{D}}}_{{conv}}$ in (3.8). In our cases, the liquid phase contribution to (3.8) is almost independent of 
$Sc_g$. The influence of 
$Sc_g$ on 
${\boldsymbol{\mathsf{D}}}_{{conv}}$ is evident in the cases with 
$\nabla ^v \langle c \rangle$ since the bubbles have an average velocity in the vertical direction due to buoyancy. The high average gas phase velocity in the gradient direction induces large average gas phase scalar fluctuations. The high gas phase velocity and scalar fluctuations cause a significant gas phase contribution 
$\phi \langle \boldsymbol {u}'c' \rangle _g$ to 
${\boldsymbol {D}}^{v}_{{conv}}$ although the parameter 
$\phi$ is small.
The proposed scalings for 
${\boldsymbol {D}}_{{conv}}$ are shown in figure 8(c) with a proportionality constant of 
$0.94$ for 
${\boldsymbol {D}}^v_{{conv}}$ and 
$0.25$ for 
${\boldsymbol {D}}^h_{{conv}}$. In the top panel, the dash-dotted line represents the scaling for 
${\boldsymbol {D}}^v_{{conv}}$ using the lowest 
$D_{{mol,s}}$ (case A4) and the dashed line is the same scaling using the highest 
$D_{{mol,s}}$ (case B1). All other cases with 
$Sc_g=0.7$ fall between these limits. Because of the previously discussed effects of 
$Sc_g$ on 
${\boldsymbol {D}}^v_{{conv}}$, the results for 
$Sc_g=7$ (grey crosses) are not in such good agreement with the proposed scalings. More data are, however, needed to study the effects of 
$Sc_g$ on 
${\boldsymbol {D}}^v_{{conv}}$. The bottom panel of figure 8(c) shows the scaling for 
${\boldsymbol {D}}^h_{{conv}}$ and all cases show that the effects of molecular diffusion (
${Sc}$) are not as significant.
4. Conclusions
We have performed DNS of bubbly flows with passive scalars and shown a transition of the scalar spectra from a 
$k^{-5/3}$ scaling (experienced in single-phase isotropic turbulence) to a 
$k^{-3}$ scaling (in the case of a mean scalar gradient in the vertical direction). In the investigated parameter ranges, we find that the transition length scale is comparable to or below the bubble diameter (
$k \geq k_{d_b}$) and decreases with the liquid scalar diffusivity (increasing 
$Sc_l$). We compute the scalar spectrum budget and show that the scalar fluctuations are produced by the mean scalar gradient at the scales above the bubble diameter (
$k < k_{d_b}$), while the production term decays as 
$k^{-3}$ for 
$k > k_{d_b}$. This is opposed to the velocity fluctuations that are continuously produced and directly dissipated at scales 
$k > k_{d_b}$.
The observed scalar production behaviour is valid in the presence of a mean scalar gradient. However, if the scalar in the liquid is injected/consumed at the bubble surfaces, we expect a scalar production peak at 
$k_{d_b}$ and hence a production term more similar to the momentum transfer. At length scales below the bubble diameter, the scalar transfer term in the budget equation shows a 
$k^{-1}$ scaling in a narrow range of wavenumbers inducing the 
$k^{-3}$ scaling of the scalar spectra for cases with a mean scalar gradient in the vertical direction. These findings are in agreement with the hypothesis proposed in Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023) about the physical mechanisms behind the 
$k^{-3}$ scaling found in that experimental work for heat transport in bubbly flows with an active turbulent grid.
We also examined the scalar effective convective diffusivity and the Sherwood number of the bubble suspension and find that 
${Sh} \propto {Pe}_\phi$ (and consequently 
${\boldsymbol {D}}^h_{{conv}} \propto \sqrt {g d_b}\phi ^{0.4}d_b$) for an imposed scalar gradient in the horizontal direction. For a vertical scalar gradient, we find, however, 
${Sh} \propto {Pe}^{1.15}_\phi$ and 
${\boldsymbol {D}}^v_{{conv}} \propto (\sqrt {g d_b}d_b)^{1.15}\phi ^{0.46}/D_{{mol,s}}^{0.15}$ for a constant gas scalar diffusivity. These scalings are based on a priori known parameters and extend the model proposed by Alméras et al. (Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015). We find that the gas scalar diffusivity significantly influences the average gas scalar fluctuations when the scalar gradient is in the vertical direction. This effect modifies 
${\boldsymbol {D}}^v_{{conv}}$ and should be considered when developing improved models. Future investigations are needed to investigate how different parameters, such as 
$Ga$ and 
$Eo$, volume fraction, molecular diffusivities of the two phases and different scalar injection mechanisms, influence the observed scalar dynamics.
Funding
The authors acknowledge support by the Swedish Research Council (Vetenskapsrådet), grant VR 2017-05031. The authors gratefully acknowledge the HPC RIVR consortium (www.hpc-rivr.si) and EuroHPC JU (eurohpc-ju.europa.eu) for funding this research by providing computing resources of the HPC system Vega at the Institute of Information Science (www.izum.si). Other computational resources have been provided by the Swedish National Infrastructure for Computing (SNIC) at NSC partially funded by the Swedish Research Council through grant agreement no. 2018-05973.
Declaration of interests
The authors report no conflict of interest.
Open access
