# Direct numerical simulations of bubble-mediated gas transfer and dissolution in quiescent and turbulent flows

Palas Kumar Farsoiya, Quentin Magdelaine, Arnaud Antkowiak, Stéphane Popinet, Luc Deike

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Palas Kumar Farsoiya, Quentin Magdelaine, Arnaud Antkowiak, Stéphane Popinet, Luc Deike. Direct numerical simulations of bubble-mediated gas transfer and dissolution in quiescent and turbulent flows. Journal of Fluid Mechanics, 2023, 954, pp.A29. 10.1017/jfm.2022.994. hal-03933249

HAL Id: hal-03933249

https://hal.science/hal-03933249

Submitted on 10 Jan 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Direct numerical simulations of bubble-mediated gas transfer and dissolution in quiescent and turbulent flows

Palas Kumar Farsoiya<fnref n="1" />, Quentin Magdelaine<fnref n="2" />, Arnaud Antkowiak<fnref n="2" />, Stéphane Popinet<fnref n="2" /> and Luc Deike<fnref n="1" />,<fnref n="3" />,†

(Received 7 May 2022; revised 18 November 2022; accepted 20 November 2022)

We perform direct numerical simulations of a gas bubble dissolving in a surrounding liquid. The bubble volume is reduced due to dissolution of the gas, with the numerical implementation of an immersed boundary method, coupling the gas diffusion and the Navier-Stokes equations. The methods are validated against planar and spherical geometries' analytical moving boundary problems, including the classic Epstein-Plesset problem. Considering a bubble rising in a quiescent liquid, we show that the mass transfer coefficient $k_L$ can be described by the classic Levich formula $k_L = (2/\sqrt{\pi})\sqrt{\mathcal{D}_l} U(t)/d(t)$, with $d(t)$ and $U(t)$ the time-varying bubble size and rise velocity, and $\mathcal{D}_l$ the gas diffusivity in the liquid. Next, we investigate the dissolution and gas transfer of a bubble in homogeneous and isotropic turbulence flow, extending Farsoiya et al. (J. Fluid Mech., vol. 920, 2021, A34). We show that with a bubble size initially within the turbulent inertial subrange, the mass transfer coefficient in turbulence $k_L$ is controlled by the smallest scales of the flow, the Kolmogorov $\eta$ and Batchelor $\eta_B$ microscales, and is independent of the bubble size. This leads to the non-dimensional transfer rate $Sh = k_L L^* / \mathcal{D}_l$ scaling as $Sh/Sc^{1/2} \propto Re^{3/4}$, where $Re$ is the macroscale Reynolds number $Re = u_{rms}L^*/v_l$, with $u_{rms}$ the velocity fluctuations, $L^*$ the integral length scale, $v_l$ the liquid viscosity, and $Sc = v_l/\mathcal{D}_l$ the Schmidt number. This scaling can be expressed in terms of the turbulence dissipation rate $\epsilon$ as $k_L \propto Sc^{-1/2}(\epsilon v_l)^{1/4}$, in agreement with the model proposed by Lamont & Scott (AIChE J., vol. 16, issue 4, 1970, pp. 513-519) and corresponding to the high $Re$ regime from Theofanous et al. (Intl J. Heat Mass Transfer, vol. 19, issue 6, 1976, pp. 613-624).

**Key words:** bubble dynamics, coupled diffusion and flow

† Email address for correspondence: ldeike@princeton.edu

<footnote n="1">¹Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="2">²Institut Jean Le Rond d'Alembert, CNRS UMR 7190, Sorbonne Université, Paris 75005, France</footnote>

<footnote n="3">³High Meadows Environmental Institute, Princeton University, Princeton, NJ 08544, USA</footnote>
# 1. Introduction

Gas transfer from bubbles to liquids is critical in multiple chemical engineering and environmental processes such as scrubbers, bubble reactors and lake oxygenation, and depends on the local flow agitation (Risso 2018), while at the ocean surface, air entrainment by breaking waves contributes to the exchange of gases critical for the climate system (Garbe et al. 2014; Reichl & Deike 2020; Deike 2022). Modelling approaches of bubble-mediated gas transfer in a turbulent environment require knowledge of the individual bubble mass transfer coefficient $k_L$ (Woolf & Thorpe 1991; Keeling 1993; Liang et al. 2011; Deike & Melville 2018).

Mass exchange related to bubble dissolution presents similarities to other heat and mass transfer problems. The classic Stefan problem (Stefan 1891) describes heat or mass transfer coupled with a moving boundary and was initially studied in the context of ice growth at solid-liquid interfaces (Crepeau 2007). Akin to the spherical Stefan problem, Epstein & Plesset (1950) provided a quasi-stationary analytical solution to the stagnant spherical bubble dissolution. Peñas-López et al. (2016) studied the bubble growth and dissolution by removing the quasi-static radius approximation, and investigated the effect of the history term that relates to the mass transfer coupled with moving boundary. Lohse & Zhang (2015) and Zhu et al. (2018) investigated a surface nano-bubble experimentally and numerically, comparing to the Epstein & Plesset (1950) results, and found that pinning of the contact line and oversaturation stabilizes the nano-bubble from dissolution.

The mass transfer from a rising bubble in a steady state (ignoring the moving boundary effect) was described by Levich (1962) using a thin boundary layer and potential flow assumptions around the bubble, leading to a mass transfer coefficient given by $k_L = (2/\sqrt{\pi})\sqrt{D_l}U/d$, with $d$ the bubble diameter, $U$ the bubble rise velocity, and $D_l$ the gas diffusivity in the liquid. This can be written in terms of the non-dimensional transfer rate (Sherwood number) $Sh = k_Ld/D_l$ as $Sh = (2/\sqrt{\pi}) Pe^{1/2}$, where $Pe = Re Sc$ is the Péclet number, $Re = Ud/v_l$ is the Reynolds number, and $Sc = v_l/D_l$ is the Schmidt number, with $v_l$ the liquid viscosity. These approximations are similar to the penetration theory (Higbie 1935) that is based on the fact that in unsteady heat or mass transfer, the depth of penetration into the exposed boundary is dependent on the time of contact. The longer the contact time, the deeper the penetration, and the equivalence for both theories is valid for high Péclet numbers (ratio of the rate of advective transport to the diffusive transport; Sideman 1966). The application of the penetration theory then relies on finding the correct time scale of contact. An experimental study on rising CO₂ bubbles by Bowman & Johnson (1962) found that the mass transfer coefficients were in close agreement with the theoretical predictions from Higbie (1935) and Levich (1962). Takemura & Yabe (1998) developed an experimental set-up to measure the size of individual O₂ bubbles rising in silicone oil, and estimated the mass transfer coefficient with results comparable to the penetration theory. The bubbles formed at the atmosphere-ocean interface correspond to high Schmidt number ($Sc \approx 600$ for CO₂ at 20 °C), while the $Sc \ll 1$ regime occurs in liquid metals and astrophysics (Gotoh et al. 2013).

A configuration of practical importance in industrial applications involves dense bubble clouds, or swarms, and requires estimations of the bubble gas transfer coefficient in this bubbly turbulent environment. Roghair (2012) found that the gas transfer coefficient increases with gas hold-up (the ratio of gas volume to total volume) by numerical simulations, while experimental studies on bubble swarms (Colombet et al. 2011, 2015) revealed that the global gas transfer coefficient at high Schmidt numbers is described
approximately by the Levich (1962) formula $k_L = (2/\sqrt{\pi})\sqrt{\mathcal{D}_l\bar{U}/d}$, considering the mean bubble rise speed $\bar{U}$ in the swarm, which is affected by turbulent fluctuations.

The mass transfer of a passive scalar at the surface of a turbulent flow has received much attention. Fortescue & Pearson (1967) studied the absorption of CO₂ in a turbulent water flow through an open channel, and approximated the mass transfer coefficient using a large eddy model $k_L \propto (\mathcal{D}_l u_{rms}/L^*)^{1/2}$, where $u_{rms}$ is the root-mean-square velocity of the turbulent flow. Theofanous, Houze & Brumfield (1976) identified two distinct mass transfer regimes for energy-containing and energy-dissipating turbulent motions, i.e. $k_L \propto (\mathcal{D}_l u_{rms}/L^*)^{1/2}$ for low Re, and $k_L \propto (v_l/\mathcal{D}_l)^{-1/2}(\epsilon v_l)^{1/4}$ at high Re, where $\epsilon$ is the turbulence dissipation rate. Here, the Reynolds number is defined using a macroscopic scale allowing for comparisons between various configurations. These scalings can be written in terms of the non-dimensional transfer rate $Sh = k_L L^*/\mathcal{D}_l$, where $L^*$ is the integral length scale, as $Sh \propto (Sc Re)^{1/2}$ (for low Re) and $Sh \propto Sc^{1/2}Re^{3/4}$ (for high Re). Herlina & Wissink (2016, 2019) and Pinelli et al. (2022) estimated the critical Reynolds number where this transition occurs as $Re \approx 500$ using numerical simulations. Lamont & Scott (1970) presented experiments on mass transfer from bubbles in concurrent turbulent flows, and argued that the small-scale motions of the turbulent flow control the mass transfer, leading to $k_L \propto (v_l/\mathcal{D}_l)^{-1/2}(\epsilon v_l)^{1/4}$, corresponding to the high Re regime from Theofanous et al. (1976). An earlier discussion of a gas bubble suspended in a turbulent liquid stream was provided by Levich (1962), who proposed an estimate for the gas transfer coefficient from a single bubble similar to the low Re regime from Theofanous et al. (1976), and mentioned the difficulty of experimental verification of the theories. In Farsoiya, Popinet & Deike (2021), we performed direct numerical simulations of a single bubble in homogeneous isotropic turbulence for a dilute gas component (so that the bubble size did not change as it transferred gas) and bubble size in the inertial subrange. We proposed that the gas transfer coefficient scaled as $Sh \propto (Sc Re)^{1/2}$, comparable to the low Re regime from Theofanous et al. (1976), with simulations performed at relatively low Re, without variations of the bubble size.

The importance of the topic has motivated the development of numerical methods to solve mass transfer in bubbly flows. Sato, Jung & Abe (2000) and Davidson & Rudman (2002) developed numerical techniques where the gas concentration is continuous across the fluid–fluid interface, while the effect of the discontinuous concentration due to solubility is considered by Bothe et al. (2004). Darmana, Deen & Kuipers (2006) provided a three-dimensional front tracking model with mass transfer. Haroun, Legendre & Raynal (2010) and Marschall et al. (2012) presented a one-fluid formulation for the algebraic volume of fluid method. Bothe & Fleckenstein (2013) introduced a two-field approach using a geometrical volume of fluid method for multicomponent conjugate mass transfer. Subgrid scale models to simulate high-Schmidt-number bubble-mediated mass exchange have been developed by Bothe & Fleckenstein (2013), Weiner & Bothe (2017) and Claassen et al. (2020). All the above numerical techniques do not consider local changes in volume due to gas dissolution. Recent advances in numerical methods by Tanguy et al. (2014), Fleckenstein & Bothe (2015), Dodd (2017), Scapin, Costa & Brandt (2020), Maes & Soulaine (2020) and Gennari, Jefferson-Loveday & Pickering (2022) allow us to simulate problems of mass transfer with local volume changes.

In the present study, the initial size of the bubble is in the inertial range of the turbulent flow, and gas from the bubble dissolves in the liquid outside as the bubble volume decreases. We use an immersed boundary method based on volume penalization (Magdelaine 2019) to model the gas diffusion from the bubble. This paper is organized as
follows. In § 2, we present the numerical framework and validate our approach against exact solutions of the moving interface due to dissolution in the planar and spherical geometries, and discuss the special case of the Epstein & Plesset (1950) problem. In § 3, we apply our numerical methods to study the mass transfer coefficient of a bubble dissolving while rising in quiescent water, and show that the classic Levich formula can describe the mass transfer coefficient, extended for a time-varying bubble size $d(t)$ and velocity $U(t)$, leading to $k_L = (2/\sqrt{\pi})\sqrt{\mathcal{D}_l}U(t)/d(t)$. In § 4, we study the dissolution of a single bubble in an otherwise homogeneous and isotropic turbulence flow, and show that the mass transfer coefficient can be described by $Sh/Sc^{1/2} \propto Re^{3/4}$, with the turbulent Sherwood number $Sh = k_L L^*/\mathcal{D}_l$ and macroscale Reynolds number $Re = u_{rms}L^*/v_l$, with $u_{rms}$ the velocity fluctuations, and $L^*$ the integral length scale. This relationship is equivalent to $k_L \propto Sc^{-1/2}(\epsilon v_l)^{1/4}$, showing that the smallest scales control the gas transfer in the flow, namely the Kolmogorov and Batchelor length scales.

## 2. Numerical framework

### 2.1. The Basilisk solver

We solve the three-dimensional (3-D), incompressible, two-phase Navier-Stokes equations using the open-source solver Basilisk (Popinet 2009, 2013, 2015):

$$ \frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot (\mathbf{u}\mathbf{u}) = \frac{1}{\rho} [-\nabla p + \nabla \cdot (\mu(\nabla \mathbf{u} + \nabla \mathbf{u}^T))] + \frac{\gamma}{\rho} \kappa \delta_s \mathbf{n}, \quad (2.1) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (2.2) $$

$$ \frac{\partial \mathcal{T}}{\partial t} + \mathbf{u} \cdot \nabla \mathcal{T} = 0, \quad (2.3) $$

where $\mathbf{u}$, $p$, $\gamma$, $\mu$, $\rho$, $\kappa$, $\mathbf{n}$ and $\mathcal{T}$ are the velocity, pressure, surface tension coefficient, viscosity, density, curvature, interface normal and volume fraction fields, respectively. The Dirac distribution function $\delta_s$ concentrates the surface tension force on the interface (Popinet 2009). The solver has been validated extensively for complex interfacial flows (Popinet 2015; van Hooft et al. 2018; Berny et al. 2020; Gumulya et al. 2020; Farsoiya et al. 2021; Mostert, Popinet & Deike 2021; Rivière et al. 2021). It uses the projection method to compute the velocity and pressure, and the geometric volume of fluid method for the evolution of the interface between two immiscible fluids (Tryggvason, Scardovelli & Zaleski 2011). The piecewise linear interface calculation geometric interface and flux reconstruction ensures a sharp representation of the interface (Scardovelli & Zaleski 1999; Marić, Kothe & Bothe 2020) and is combined with an accurate height-function curvature calculation and a well-balanced, continuum surface tension model (Brackbill, Kothe & Zemach 1992; Popinet 2018). Well-balanced numerical schemes ensure that specific equilibrium solutions of the (continuous) partial differential equations are recovered by the (discrete) numerical scheme. Such schemes are important when surface tension is the dominant force, such as in two-phase interfacial flows including drops or bubbles close to Laplace's equilibrium (Popinet 2018). Test cases by Popinet (2009) demonstrate second-order convergence towards the exact circular equilibrium shape as a function of spatial resolution.
## 2.2. Immersed boundary method for gas diffusion

We use the numerical method developed by Magdelaine (2019), which relies on a distributed Dirac-like forcing term, as in the immersed boundary method of Peskin (1972), in order to enforce Dirichlet conditions for the gas concentration on the evolving interface. The time evolution of the gas concentration is given by

$$ \frac{\partial c}{\partial t} + \nabla \cdot (\mathbf{u}c) = \nabla \cdot (\mathcal{D}_l \nabla c) + F, \quad (2.4) $$

where the forcing term $F$ relaxes the concentration toward its gas phase value on a time scale $\tau_c$ in the gas phase region $B$, and can be written as a surface integral:

$$ F(x) = \int_{x(s) \in B} \left( \frac{c_g - c}{\tau_c} \right) \delta(x - x(s)) \, ds, \quad (2.5) $$

where $\delta(x)$ is a two-dimensional Dirac function, and $ds$ represents the surface element of the boundary (Peskin 1972). Equation (2.4) can then be re-expressed as

$$ \frac{\partial c}{\partial t} + \nabla \cdot (\mathbf{u}c) = \nabla \cdot (\mathcal{D}_l \nabla c) - \frac{c - c_g}{\tau_c} \mathcal{H}(x). \quad (2.6) $$

Here, $\mathcal{H}(x) = \int_{-\infty}^{0} \delta(x - x(s)) \, ds$ is a (surface) Heaviside function for the non-diffusive phase (gas phase) with the relaxation time scale $\tau_c = 10^{-6} \Delta^2 / \mathcal{D}_l$, with $\Delta$ the grid size. The relaxation time scale needs to be small enough to ensure an accurate imposition of the Dirichlet condition on the interface, while being large enough to avoid possible amplification of round-off errors.

The assumption of continuous chemical potentials, which is good for most applications at an interface $\Sigma(t)$, results in Henry's law (Bothe & Fleckenstein 2013):

$$ c_l = c_g \alpha \quad \text{at} \quad \Sigma(t), \quad (2.7) $$

where the dimensionless ratio of the liquid phase concentration to the gas phase concentration of the component transferred, $\alpha$, is Henry's law solubility constant (solubility hereafter). We consider a single component and an isothermal configuration, and there is no bubble break-up or large deformation of the interface that can change the pressure inside the bubble significantly. Hence the solubility, which depends on temperature, pressure and mole fraction of the gas components (Bothe & Fleckenstein 2013), is assumed to be constant. We modify the last term in (2.4) to account for the solubility:

$$ \frac{\partial c}{\partial t} + \nabla \cdot (\mathbf{u}c) = \nabla \cdot (\mathcal{D}_l \nabla c) - \frac{c - \mathcal{E}}{\tau_c} \mathcal{H}, \quad (2.8) $$

where $\mathcal{E}$ is a conditional variable equal to $\alpha c_g$ and $c_g$ in the interfacial and non-interfacial cells, respectively. The gradients of $c$ are calculated using central differences, and the Heaviside function is approximated as $1 - \mathcal{T}$.

The advection and diffusion terms in (2.8) are treated in a time-split manner. For diffusion we use a time-implicit Euler discretization

$$ \frac{c^{n+1} - c^n}{\Delta t} = \nabla \cdot (\mathcal{D}_l \nabla c^{n+1}) + \beta c^{n+1} + r, \quad (2.9) $$

where $\beta = -\mathcal{H}/\tau_c$ and $r = -\mathcal{H}\mathcal{E}/\tau_c$. Equation (2.9) gives a set of linear equations that are solved efficiently using the multigrid method (Popinet 2015).
The solubility boundary condition (2.7) presents a discontinuity for the concentration field across the interface similar to the volume fraction field $T$. As discussed by Bothe & Fleckenstein (2013), when the mass transfer is coupled with fluid flow, the concentration discontinuity given by Henry's law must be satisfied at the fluid-fluid interface. Moreover, the advection of the concentration field $c$ must be consistent with the advection of the volume fraction field since any relative motion would lead to artificial mass transfer across the interface. To avoid this numerical diffusion across the interface, we use the two-scalar approach for advection proposed by López-Herrera et al. (2015), i.e. two tracer fields $\phi_g$ and $\phi_l$ associated with the volume of fluid $T$ are defined:

$$ \phi_l = c_l T, \quad \phi_g = c_g(1-T), \qquad (2.10a,b) $$

$$ c = Tc_l + (1-T)c_g. \qquad (2.11) $$

Using (2.7) and (2.11), we get

$$ \phi_l^n = \frac{c^n \alpha T}{\alpha T + (1-T)}, \quad \phi_g^n = \frac{c^n (1-T)}{\alpha T + (1-T)}. \qquad (2.12a,b) $$

The advection equation for $\phi_{g/l}$ reads

$$ \frac{\partial \phi_{g/l}}{\partial t} + \nabla \cdot (u \phi_{g/l}) = 0, \qquad (2.13) $$

and is solved using the volume-of-fluid associated fields (López-Herrera et al. 2015; Farsoiya et al. 2021) that guarantee strictly non-diffusive transport close to the interface. The concentration tracer is updated after advection using

$$ c = \phi_g + \phi_l. \qquad (2.14) $$

We then calculate the dissolution velocity at the interface, $u^\Sigma$, which is used to update the volume fractions (Magdelaine 2019):

$$ \frac{\partial T}{\partial t} + u^\Sigma \cdot \nabla T = 0. \qquad (2.15) $$

The boundary condition for molar mass flux at an interface is given by (Fleckenstein & Bothe 2015)

$$ [[c(u - u^\Sigma)]] \cdot n + [[j]] \cdot n_\Sigma = 0, \qquad (2.16) $$

where $u$ and $u^\Sigma$ are the barycentric and interface velocities, respectively.

The Stefan currents from a dissolving bubble can be estimated from the source term in the continuity equation. They are discussed by Dodd (2017) in the context of evaporating droplets in homogeneous and isotropic turbulence (HIT) flow

$$ \nabla \cdot u = \frac{1}{Pe} \left( \frac{\nabla Y_v \cdot n \delta(n)}{1 - Y_v} \right) \left( 1 - \frac{\rho_g}{\rho_l} \right), \qquad (2.17) $$

where $Y_v$ is the mass fraction of the vapour. The numerical simulations performed for the bubble in turbulent flow are in the region $\rho_g \ll \rho_l$ ($\rho_l/\rho_g \approx 850$) and $Pe \gg 1$ ($80 \lesssim Pe \lesssim 5 \times 10^4$), the Péclet number being defined based on some characteristic velocity of the flow. Hence the Stefan currents are small compared to the flow around the bubble. In addition, the vapour recoil and surface tension typical of air bubbles dissolution in water scale as $(\rho_g(u - u^\Sigma) \cdot n_\Sigma)^2(1/\rho_g) \approx O(10^{-6})$ and $\gamma\kappa \approx O(10^2)$, respectively.
These forces can be comparable for the cases of strong evaporation of a droplet and cavitation of a bubble (Fleckenstein & Bothe 2015), which is not the case in our study. Consequently, Stefan currents and vapour recoil are not included in the present numerical framework.

As the Stefan currents are not modelled, the interface velocity $u^\Sigma$ from (2.16) reduces to

$$u^\Sigma = \frac{M_w}{\rho_g(1-\alpha)} \mathcal{D}_l \nabla c, \qquad (2.18)$$

where $M_w$ is the molecular mass of the gas. The time step is not restricted by the scalar diffusion terms as we are using an implicit discretization for diffusion terms. However, fast diffusion calculates the phase change interface velocity, which restricts the time step, and the interfacial velocity $u^\Sigma$ restricts the computational time step in accordance with Courant–Friedrichs–Lewy criteria.

Moreover, this approach is convenient in simulating HIT flows as these require periodic conditions on the computational domain boundaries, and having a source term in the continuity equation is then needed to be compensated for the extra liquid inside the domain. We thus assume that the solvent is an infinite reservoir, and study the bubble dissolution in an HIT flow. Note also that there are methods developed by Fleckenstein & Bothe (2015), Malan (2018), Scapin et al. (2020) and Maes & Soulaine (2020) for inclusion of Stefan currents. A recent independent development by Gennari et al. (2022) in Basilisk based on the two-scalar approach by Fleckenstein & Bothe (2015) treats the velocity discontinuity at the interface that is due to the mass transfer between the phases.

### 2.3. Validation

#### 2.3.1. Stefan problem

We test our methods with the classical moving boundary problem, called the Stefan problem (Stefan 1891) in planar geometry, where the liquid and gas phases are separated by a moving interface $Y(t)$. The concentration of gas in the liquid phase (neglecting the advection) is given by

$$\frac{\partial c}{\partial t} = \mathcal{D}_l \frac{\partial^2 c}{\partial y^2}, \quad \text{for } y > Y(t) \quad (\text{liquid phase}). \qquad (2.19)$$

The moving interface concentration is constant:

$$c(Y(t), t) = \alpha c_g, \quad t \ge 0, \qquad (2.20)$$

and

$$c(y \to \infty, t) = 0, \quad c(y > Y(t), 0) = 0, \qquad (2.21a,b)$$

where $\alpha$ is Henry's solubility constant, and $c_g$ is the gas concentration in the gas phase. The problem can be non-dimensionalized using the characteristic length $L$, $\zeta = (y - Y(0))/L$, the time as a mass Fourier number $Fo_m = \mathcal{D}_l/L^2t$, the moving boundary
$\Sigma(Fo_m) = (Y(t) - Y(0))/L$, and the concentration $\Theta = c(y, t)/(\alpha c_g)$, so that

$$ \frac{\partial \Theta}{\partial Fo_m} = \frac{\partial^2 \Theta}{\partial \zeta^2}, \quad \text{for } \zeta > \Sigma(Fo_m) \quad \text{(liquid phase).} \qquad (2.22) $$

The moving interface concentration is constant:

$$ \Theta(\Sigma, Fo_m) = 1, \quad Fo_m \ge 0. \qquad (2.23) $$

The velocity of the moving interface from (2.16), $Y(t)$, is given by

$$ \frac{d\Sigma}{dFo} = St \left(\frac{\partial\Theta}{\partial\zeta}\right)_{\zeta=\Sigma(Fo_m)}, \quad Fo_m > 0, \qquad (2.24) $$

with conditions

$$ \Theta(\zeta \to \infty, t) = 0, \quad \Sigma(0) = 0, \quad \Theta(\zeta > 0, 0) = 0, \qquad (2.25a-c) $$

and Stefan number

$$ St = \frac{\alpha}{1 - \alpha}. \qquad (2.26) $$

The time evolution of the interface and concentration is given by (Alexiades & Solomon 2018)

$$ \Sigma(Fo_m) = 2\lambda\sqrt{Fo_m}, \qquad (2.27) $$

$$ \Theta(\zeta, Fo_m) = \frac{1}{\operatorname{erfc}(\lambda)} \operatorname{erfc}\left(\frac{\zeta}{2\sqrt{Fo_m}}\right), \qquad (2.28) $$

where $\lambda$ is the solution of the equation

$$ \lambda \operatorname{erfc}(\lambda) e^{\lambda^2} = -\frac{St}{\sqrt{\pi}}. \qquad (2.29) $$

Comparing the integration of concentration profile (2.28),

$$ \int_{\Sigma(Fo_m)}^{\infty} \Theta(\zeta, Fo_m) d\zeta = \int_{0}^{\infty} \left[ \frac{1}{\operatorname{erfc}(\lambda)} \operatorname{erfc}\left(\frac{\zeta}{2\sqrt{Fo_m}}\right) \right] d\zeta = -2\lambda\sqrt{Fo_m}, \qquad (2.30) $$

with (2.27) shows that the solution conserves mass. Note that (2.27) (displacement of interface) and (2.30) (mass of gas) have negative and positive signs, respectively ($\lambda < 0$).

The test simulation is performed in a two-dimensional square box with $St = 0.8$, with an adaptive resolution of 2048 cells per side of the domain. The position of the interface with time is shown in figure 1(a), and the concentration of gas in the liquid at a location in the domain is shown in figure 1(b). Both quantities show good agreement with the theoretical predictions (2.27) and (2.28), respectively. In figure 1(c), the net mass of gas is conserved, and liquid side mass shows good agreement with (2.30). Figures 1(d,e) show first-order convergence of the numerical methods. The agreement is good throughout the full time and space evolution.
<figure><img src="image_4.png" /><figcaption>Bubble-mediated gas transfer</figcaption></figure>

Figure 1. Test results for the planar Stefan problem, for $St = 0.8$ and grid size $(2^{11})^2$ (adaptive). (a) The interface position shows good agreement with (2.27). (b) Concentration at $\zeta = 0.02$ shows good agreement with (2.28). (c) Gas mass is conserved and shows good agreement with (2.30). The scripts for the test cases are available; see Farsoiya, Popinet & Deike (2022a). (d) Maximum relative error at different resolutions, $\max|c_{11} - c_N|/c_s$, where $c_{11}$ and $c_N$ are numerical solutions at resolution $2^{11}$ and lower, respectively, displaying first-order convergence. (e) Maximum relative error in mass conservation at different resolutions, $\max|m(t) - m(0)|/m(0)$, where $m(t)$ and $m(0)$ are net mass (time evolution) and initial net mass, displaying first-order convergence.

### 2.3.2. Dissolution of a static spherical bubble

We now apply our numerical methods to the dissolution of a static spherical bubble. The gas diffusion in the liquid phase outside a spherical bubble, neglecting advection,
is governed by

$$ \frac{\partial c}{\partial t} = \mathcal{D}_l \frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 \frac{\partial c}{\partial r} \right), \quad \text{for } r > R(t) \quad \text{(liquid phase).} \qquad (2.31) $$

The moving interface concentration is constant:

$$ c(R(t), t) = \alpha c_g, \quad t \ge 0, \qquad (2.32) $$

and $c(r \to \infty, t) \to 0$, where $\alpha$ is Henry's solubility constant, and $c_g$ is the gas concentration in the gas phase. The velocity of the moving interface $R(t)$ from (2.16) is given by

$$ \frac{dR(t)}{dt} = \frac{\mathcal{D}_l}{c_g(1-\alpha)} \left(\frac{\partial c}{\partial r}\right)_{r=R(t)}, \quad t > 0, \qquad (2.33) $$

with initial conditions $R(0) = R_0$ and $c(r > R_0, 0) = 0$. We consider the configuration for a large driving force for transport given by a large Stefan number (which will apply in particular to larger solubility), here $St = 0.8$. At this Stefan number, the change in bubble size affects the diffusion of gas ($St \sim O(1)$) over time. We solve numerically (forward Euler method) the moving boundary problem given by (2.31)–(2.33) (numerical moving boundary problem, NMBP). The concentration $c$ and bubble radius $R(t)$ can be non-dimensionalized as $\Theta = c(r, t)/(\alpha c_g)$ and $\epsilon = R(t)/R_0$ (2.36a-c).

The non-dimensional bubble radius is shown in figure 2(a), and the concentration at a location outside the bubble at $r = R_0 + 0.2$ is shown in figure 2(b). Good agreement is found with the numerical solution to the moving boundary problem. Figure 2(c) illustrates the mass conservation properties. The numerical method is robust in keeping the shape of the shrinking bubble nearly spherical, as shown in figure 2(d), with the sphericity $\Psi = \pi^{1/3}(6V_b)^{2/3}/A_b$, where $V_b$ and $A_b$ are the volume and area of the bubble, respectively. As the bubble shrinks, the number of cells per diameter decreases, and the amplitude of the oscillations in sphericity increases, which is expected.

### 2.3.3. Epstein-Plesset test

Epstein & Plesset (1950) provided an analytical solution on the dissolution of a spherical bubble under the quasi-static radius approximation (Weinberg & Subramanian 1980; Peñas-López et al. 2016), i.e. treating the bubble radius as a constant at the concentration boundary condition at the interface. It is valid if the motion of the interface is much slower than the diffusion process. The Epstein & Plesset (1950) test and approximations have been discussed in Duncan & Needham (2004), Peñas-López et al. (2016), Chu & Prosperetti (2016) and Zhu et al. (2018).

Under these approximations, the time evolution of the interface is given by Epstein & Plesset (1950) as

$$ \frac{dR}{dt} = -\mathcal{D}_l \frac{M_w(\alpha c_g - c_0)}{\rho_g} \left\{ \frac{1}{R} + \frac{1}{\sqrt{\pi \mathcal{D}_l t}} \right\}. \qquad (2.34) $$

In dimensionless form,

$$ \frac{d\epsilon}{dx} = -\frac{x}{\epsilon} - 2\xi, \qquad (2.35) $$
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 2. Dissolution of a static bubble in the high-St-regime (St = 4) large driving force for transport. Grid size 200 cells per initial diameter. (a) Time evolution of bubble radius. (b) Concentration at r = R₀ + 0.2 shows good agreement with numerical solution of the moving boundary problem given by (2.31)–(2.33). (c) The net mass of gas is conserved, and the evolution of the mass of gas dissolved in the liquid agrees with the volume integration of c (NMBP). (d) Sphericity remains close to unity as the bubble shrinks.</figcaption></figure>

where

$$ \epsilon = \frac{R(t)}{R_0}, \quad x^2 = \frac{2\mathcal{D}|M_w(\alpha c_g - c_0)}{\rho_g R_0^2} t, \quad \xi = \sqrt{\frac{M_w(\alpha c_g - c_0)}{2\pi\rho_g}} = \sqrt{\frac{St}{2\pi}}, \quad (2.36a-c) $$

where St is the Stefan number. The solution of (2.35) is given in the parametric form

$$ \epsilon = e^{-\xi z} \left( \cos \left( z \sqrt{1 - \xi^2} \right) - \xi \frac{1}{\sqrt{1 - \xi^2}} \sin \left( z \sqrt{1 - \xi^2} \right) \right), \quad (2.37a) $$

$$ x = \frac{1}{\sqrt{1 - \xi^2}} e^{-\xi z} \sin \left( z \sqrt{1 - \xi^2} \right). \quad (2.37b) $$

The analytical solution of (2.31) (Darmana et al. 2006) can be used for a quasi-static boundary R(t) to get the concentration as

$$ c = c_0 + (\alpha c_g - c_0) \frac{R(t)}{r} \text{erfc} \left( \frac{r - R(t)}{\sqrt{4\mathcal{D}|t|}} \right). \quad (2.38) $$

We set up a test case to mimic the Epstein & Plesset (1950) conditions by considering $u_{EP}^\Sigma = 2 \times 10^{-4} u^\Sigma$ in (2.18), valid for a quasi-stationary bubble. This approximation
<figure><img src="image_9.png" /><figcaption>Figure 3. Dissolution of a static bubble in the quasi-stationary approximation $u_{EP}^{\Sigma} = 2 \times 10^{-4}u^{\Sigma}$. The grid size is 200 cells per initial diameter. (a) The time evolution of the bubble radius shows good agreement with (2.37a), (2.37b). (b) Concentration $\Theta = (c(r, t) - \alpha c_g)/(c_0 - \alpha c_g)$ at $r = R_0 + 0.2$, in good agreement with (2.38). Scripts for test cases are available; see Farsoiya, Popinet & Deike (2022b).</figcaption></figure>

creates a mass source and hence does not conserve the total mass of the gas. The test is performed with an adaptive local resolution of 200 cells per initial diameter. The diffusivity $\mathcal{D}_l$ non-dimensionalizes the time as a Fourier number $F_{om} = \mathcal{D}_l/R_0^2 t$, while the non-dimensional spatial variable can be written as $x^2 = 2 F_{om} St$. The motion of the bubble boundary in figure 3(a) shows good agreement with (2.37a), (2.37b). Figure 3(b) shows good agreement of the concentration at $r = R_0 + 0.2$ with (2.38). The good agreement is observed for most of the time evolution of the bubble dissolution, with small deviations at later times related to the quasi-stationary approximation by Epstein & Plesset (1950). This validates our numerical implementation in the configuration considered by Epstein & Plesset (1950).

### 3. Dissolution of a rising bubble

After considering the cases of static bubbles dissolving in a liquid, we now consider the effects of advection and the case of a bubble rising and dissolving in an initially quiescent liquid. The bubble is initially spherical and rises due to buoyancy while shrinking due to the dissolution of gas. Levich (1962) proposed a model for the gas transfer velocity from a bubble rising with terminal velocity $U$ as $k_L = 2\sqrt{\mathcal{D}_l U/\pi d}$, assuming a spherical shape and a constant size and surface concentration. Legendre & Magnaudet (1999) obtained a similar result for a steady axisymmetric straining flow mass transfer coefficient $Sh = 2\sqrt{Pe_0/\pi}$. This relationship considers a steady state with a constant terminal velocity and does not account for changes in bubble size; hence it is valid in the limit of exchange of dilute components from the bubble, and we have verified this result previously for a constant bubble size (Farsoiya et al. 2021).

We propose to extend the formula from Levich (1962) to the case of a bubble changing its size and associated rising velocity, writing

$$k_L(t) = \frac{2}{\sqrt{\pi}} \sqrt{\frac{\mathcal{D}_l U(t)}{d(t)}}, \quad (3.1)$$

where $U(t)$ and $d(t)$ are the now time-varying rise velocity and equivalent bubble diameter of a sphere having the same volume.
We compare this model against numerical simulations of a bubble rising in an initially quiescent liquid and shrinking in size due to dissolution, leading to a time-varying rise velocity $U(t)$, and instantaneous size $d(t)$. We define the non-dimensional mass transfer coefficient (or Sherwood number) $Sh = k_L(t)d_0/\mathcal{D}_l$, while the specific bubble rise configuration can be defined by the bubble Morton number $Mo = g\mu_l^4/(\rho_l\gamma^3)$ and Bond number $Bo = \rho_l g d_0^2 \gamma$ (Moore 1965; Maxworthy et al. 1996; Clift, Grace & Weber 2005; Cano-Lozano et al. 2016).

Figures 4(a–d) show a bubble rising for given Morton and Bond numbers. The bubble is resolved with up to 200 cells per diameter. The values of $Mo$ and $Bo$ are chosen from Farsoiya et al. (2021) and Deising, Marschall & Bothe (2016). In figures 4(a–d), we observe the bubble rising and shrinking in size, as well as the associated wake of gas left behind.

The total mass $N$ and average gas concentration $\bar{c}$ inside and outside the bubble is computed as

$$N_g = \int_{V_g} c \, dV, \quad N_l = \int_{V_l} c \, dV, \quad \bar{c}_g = \frac{1}{V_g} \int_{V_g} c \, dV, \quad \bar{c}_l = \frac{1}{V_l} \int_{V_l} c \, dV, \quad (3.2a-d)$$

where $V_g$ and $V_l$ are the volumes of bubble and liquid, respectively. Given that the gas does not leave the domain, the mass transfer coefficient $k_L$ is calculated as

$$k_L = \frac{1}{A_g(\alpha c_g - c_l)} \frac{dN}{dt}, \quad \text{or} \quad k_L = \frac{N_l^{n+1} - N_l^n}{A_g \Delta t (\alpha \bar{c}_g^{n+1/2} - \bar{c}_l^{n+1/2})}, \quad (3.3a,b)$$

where $A_g$ is the instantaneous surface area of the bubble.

Figures 4(e,f) show the time evolution of the non-dimensional transfer coefficient $Sh$ from the simulations. Good agreement with (3.1) is observed. Figure 4(g) shows the time evolution of the diameter of the bubble for both axisymmetric and 3-D simulations. Figure 4(h) shows that the net mass of gas is conserved as the bubble gets dissolved completely in the liquid.

We note that the gas bubble has shrunk more for the initial volume at a lower Schmidt number as it corresponds to a gas with higher diffusivity. The formula (3.1) works when provided with the instantaneous rise velocity and size of the bubble. Therefore, we expect that it would apply to a wide range of initial bubble configurations, defined by the bubble Bond and Morton numbers, including configurations at higher Reynolds numbers with spiralling or zigzagging bubble trajectories. The increase in diffusivity $\mathcal{D}_l$ increases the diffusion of gas from the bubble. Higher solubility $\alpha$ gives higher concentration at the bubble boundary and results in a higher concentration gradient and higher gas transfer. The Stefan number is $O(1)$ and does not affect (3.1) for the high Péclet regime.

## 4. Bubble dissolution in HIT flow

We now present direct numerical simulations of bubble dissolution in HIT. First, we perform a precursor simulation to achieve a stationary HIT state, and then we insert the bubble. The methodology is similar to the preceding work on bubble deformation (Perrard et al. 2021), break-up (Rivière et al. 2021) and dilute gas transfer (Farsoiya et al. 2021). A similar approach has been taken in independent studies by Loisy & Naso (2017), Dodd & Ferrante (2016) and Dodd et al. (2021). Naso & Prosperetti (2010) have observed
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

Figure 4. (a–d) Interface of a 3-D bubble (white, only left half shown in order to see the concentration field of gas inside the bubble) for various Mo and Bo (specified in the keys). (e) Non-dimensional mass transfer coefficient $Sh = k_L(t)d_0/D_l$ as a function of time $t\sqrt{g/d_0}$ for the four cases shown in (a–d). (f) Coefficient $Sh$ as a function of time for increasing $Sc = 1, 10, 100$, at $Mo = 5 \times 10^{-7}$ and $Bo = 3.125$, and corresponding curves for (3.1), using the instantaneous bubble size and velocity. (g) Time evolution of the bubble diameter for varying $Sc = 1, 10, 100$, for $Mo = 5 \times 10^{-7}$ and $Bo = 3.125$. Good agreement between the axisymmetric and 3-D simulations is observed. (h) Net gas mass is conserved until the bubble dissolves completely ($Sc = 1$, $Mo = 5 \times 10^{-7}$ and $Bo = 3.125$).
## Bubble-mediated gas transfer

an unstable behaviour in such a configuration when coupling to a rigid sphere. Such instability is specific to the two-way coupling with rigid (non-deformable) particles with a density close to that of the carrying fluid, and is not observed in deformable bubbles and droplets. The methodology is summarized briefly for completeness.

### 4.1. Precursor simulation

A cubic domain of size $L$ periodic in all three dimensions is subjected to linear volumetric forcing $f = A u(x, t)$ (Lundgren 2003) in (2.1). Rosales & Meneveau (2005) have demonstrated that this leads to a turbulent field with similar statistics to that obtained with a spectral code, and leads to a well-characterized HIT flow. We use adaptive mesh refinement on the velocity field, and the maximum level of refinement can be used to compare the resolution with that of a fixed grid. The domain consists of fluid of density $\rho_l$ and viscosity $\mu_l$. The turbulent flow is generated for increasing resolutions, with the maximum level of refinement going from 6 to 8, corresponding to an equivalent grid size of 64³ to 256³. The turbulence state is characterized by the kinetic energy density $K$, turbulence dissipation rate $\epsilon$, and Taylor microscale Reynolds number $Re_\lambda$, which are given by (Pope 2001)

$$ K = \frac{1}{V} \int_V \frac{1}{2} \rho_l |u'(x, t)|^2 dV, \quad \epsilon = \frac{1}{V} \int_V v_l \left( \frac{\partial u_i}{\partial x_j} \frac{\partial u_i}{\partial x_j} \right) dV, \quad Re_\lambda = \frac{2K}{3v_l} \sqrt{\frac{15v_l}{\epsilon}}, \quad (4.1a-c) $$

and are computed over time to characterize the turbulent flow. The root-mean-square of the velocity is $u_{rms} = \sqrt{2K/3\rho_l}$, and the eddy turnover time at the integral length scale $L^* = u_{rms}^3/\epsilon$ is given by $t_c = (L^*)^{2/3}\epsilon^{-1/3}$ (Pope 2001).

Figure 5(a) shows the evolution of the turbulent kinetic energy with time for increasing Reynolds numbers. It shows that at $t/t_c \approx 10$, the flow has reached a statistically stationary state. We show that the state is grid-independent when using an adaptive mesh refinement $L_7 \equiv 2^7$, $L_8 \equiv 2^8$ and $L_9 \equiv 2^9$. The turbulence viscous dissipation and the turbulence Reynolds number are shown in figures 5(b) and 5(c), respectively. The range of turbulence Taylor Reynolds numbers is $Re_\lambda \approx 38-150$, typical of current two-phase simulations of turbulent flow (Loisy, Naso & Spelt 2017; Elghobashi 2019; Farsoiya et al. 2021). We characterize the turbulent stationary state using the second-order structure functions in the longitudinal $D_{LL}(r)$ and transverse $D_{NN}(r)$ directions, given by $D_{LL}(r) = \frac{1}{3} \sum_i \langle (u_i(r, t) - u_i(r + d\hat{r}_i, t))^2 \rangle$ and $D_{NN}(r) = \frac{1}{6} \sum_{i \neq j} \langle (u_i(r, t) - u_i(r + d\hat{r}_j, t))^2 \rangle$, where $\hat{r}_i$ is the unit vector along the $i$th direction. Figure 5(d) shows that the scaled structure functions plateau at $C=2$ (Pope 2001) in the inertial range. The relation $D_{LL} = 3/4D_{NN}$ is verified. The bubble is inserted once the turbulent stationary state is reached, and is of a size within the inertial range. We have verified in Farsoiya et al. (2021) that the Kolmogorov length scale is reasonably well resolved in these simulations.

### 4.2. Bubble insertion

A bubble is inserted at the centroid of the cubic domain that has achieved a stationary HIT state, i.e. for $t_0 > 15t_c$. The bubble has diameter $d_0$, viscosity $\mu_g$, density $\rho_g$, and surface tension coefficient $\gamma$. The Weber number $We = 2\rho_l\epsilon^{2/3}d_0^{5/3}/\gamma \approx 2$ is kept below
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>Figure 5. Precursor simulation demonstrating the HIT flow. (a) Kinetic energy as a function of time. (b) Turbulence dissipation rate as a function of time. (c) Reynolds number at the Taylor length scale as a function of time. In all cases, these quantities reach a statistically stationary value after some time, and the bubble is inserted in this HIT flow. (d) Second-order structure function, which demonstrates good agreement with turbulence scaling.</figcaption></figure>

the critical value for break-up, $We_c$ (Rivière et al. 2021) so that the bubble might deform but never breaks. With bubble dissolution, the bubble’s diameter, and hence the Weber number, reduce further, so that we never observe bubble break-up. Note that the forcing of the turbulence continues to be applied as in Perrard et al. (2021), Rivière et al. (2021) and Farsoiya et al. (2021) once the bubble is inserted, and that we do not consider the effects of gravity/buyoancy.

The bubble has initial gas concentration $c = \rho_g/M_w$, with the solubility and diffusivity of the gas in the outside fluid given by $\alpha$ and $D_l$, respectively. The turbulence properties and simulation parameters are given in table 1. The highest resolutions for $Sc = 1$ and $Sc = 10$, 30 cases are L10 and L12, respectively. Figure 6 shows a 3-D rendering of the concentration field for $Re_\lambda = 150$, $Sc = 10$, for the time of bubble insertion (figure 6a) and then three instants during the first eddy turnover time. The bubble size can be seen as being reduced while the gas diffuses out, following turbulent fluctuations. The vorticity field is also indicated.

Another illustration of the dissolution dynamics in turbulence is provided in figure 7. Figures 7(a–d) show the reduction in bubble volume at different time instants for the case $Re_\lambda = 150$, $Sc = 1$. Note that the bubble volume reduces at a faster rate compared to the radius. The interface shown at $(t-t_0)/t_c = 2.5$ corresponding to $R(t)/R(0) \approx 0.4$ implies $V(t)/V(0) \approx 0.064$. The bubble interface is shown at an instant $(t-t_0)/t_c = 1$ in
<table><caption>Table 1. Parameters (with adaptive mesh refinement) of the simulations of bubble dissolution in turbulence. Four Reynolds numbers are used, while the Weber number (below break-up threshold), density and viscosity ratio are kept constant. The solubility $α = 0.3$ is for cases with $Sc = 1$ only.</caption><thead><tr><th>Resolution (2<sup>L</sup>)</th><th>d<sub>0</sub>/λ</th><th>L*/d<sub>0</sub></th><th>We</th><th>ρ<sub>r</sub></th><th>μ<sub>r</sub></th><th>Sc</th><th>α</th><th>Re<sub>λ</sub></th><th>Re</th></tr></thead><tbody><tr><td>L9 ≡ 2<sup>9</sup> to L12 ≡ 2<sup>12</sup></td><td>1.56</td><td>1.51</td><td>1.8</td><td>850</td><td>25</td><td>(1, 10, 30)</td><td>(0.5, 0.3)</td><td>38</td><td>85</td></tr><tr><td>L9 ≡ 2<sup>9</sup> to L12 ≡ 2<sup>12</sup></td><td>2.19</td><td>1.75</td><td>1.8</td><td>850</td><td>25</td><td>(1, 10, 30)</td><td>(0.5, 0.3)</td><td>55</td><td>170</td></tr><tr><td>L9 ≡ 2<sup>9</sup> to L12 ≡ 2<sup>12</sup></td><td>2.72</td><td>1.88</td><td>1.8</td><td>850</td><td>25</td><td>(1, 10, 30)</td><td>(0.5, 0.3)</td><td>77</td><td>400</td></tr><tr><td>L9 ≡ 2<sup>9</sup> to L12 ≡ 2<sup>12</sup></td><td>3.18</td><td>2.51</td><td>1.8</td><td>850</td><td>25</td><td>(1, 10, 30)</td><td>(0.5, 0.3)</td><td>100</td><td>800</td></tr><tr><td>L9 ≡ 2<sup>9</sup> to L12 ≡ 2<sup>12</sup></td><td>3.24</td><td>3.29</td><td>1.8</td><td>850</td><td>25</td><td>(1, 10, 30)</td><td>(0.5, 0.3)</td><td>150</td><td>1600</td></tr></tbody></table>

<figure><img src="image_22.png" /><figcaption>Figure 6. Concentration field for $Re_λ = 150$, $Sc = 10$ at (a) $t = t_0$ (bubble insertion), (b) $(t - t_0)/t_c = 0.5$, (c) $(t - t_0)/t_c = 0.75$, and (d) $(t - t_0)/t_c = 1$. Diffusion of gas around the bubble is shown by field c: the darker blue region is a high concentration inside the bubble, and yellow is a low concentration distributed due to advection and diffusion. Magnitude of vorticity is shown at the boundaries of the box behind the bubble.</figcaption></figure>

figures 7(e–h), for two $Sc$ values (1 and 10), for the same turbulence condition. The vorticity field on a plane intersecting the bubble is shown in figures 7(e, f) and the corresponding concentration field is shown in figures 7(g, h) for $Sc = 1$ and $Sc = 10$, respectively. In figure 7(h), the transient boundary layers for $Sc = 10$ are thinner compared to $Sc = 1$ (figure 7g). The smaller bubble size can also be noticed in figure 7(g) ($Sc = 1$) compared to figure 7(h) ($Sc = 10$), as more mass has diffused into the outer fluid. The concentration fields are similar to the patterns of the surrounding vorticity field.
<figure><img src="image_23.png" /><figcaption>Figure 7. Interface of bubble at different time instants showing dissolution for the case $Re_λ = 150$, $Sc = 1$, with vorticity field shown on a planar slice at the boundary behind the bubble: (a) $(t-t_0)/t_c = 0$, (b) $(t-t_0)/t_c = 1$, (c) $(t-t_0)/t_c = 2$, and (d) $(t-t_0)/t_c = 2.5$. Vorticity and concentration field (shown on a planar slice intersecting the bubble) for $Re_λ = 150$ at $(t-t_0)/t_c = 1$ for (e,g) $Sc = 1$, (f,h) $Sc = 10$.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure 8. Bubble size and non-dimensional mass transfer coefficient (based on the initial bubble size), as a function of time: (a,b) $Re_λ = 38$, $Sc = 1$; (c,d) $Re_λ = 55$, $Sc = 10$; (e,f) $Re_λ = 150$, $Sc = 10$. The $Sc = 1$ results are converged between L9 and L10, and the $Sc = 10$ results are converged between L10 and L12. Solubility $α$ rescales the constant concentration at the bubble boundary, $k_L = (dc/dt)/(αc_g)$.</figcaption></figure>

### 4.3. Dynamics and mass transfer associated with bubble dissolution in turbulence

The time evolution of the bubble radius and associated non-dimensional mass transfer coefficient are shown in figure 8 for increasing Reynolds number ($Re_λ = 38$, $Sc = 1$, $α = 0.5$, 0.3 in figures 8(a,b), $Re_λ = 55$, $Sc = 10$, $α = 0.5$ in figures 8(c,d), and $Re_λ = 150$, $Sc = 10$, $α = 0.5$ in figures 8(e,f)). The mass transfer coefficient $k_L$ is calculated using (3.3b), and the associated instantaneous $Sh(t)$ is defined as $Sh(t) = k_L(t)L^*/D_l$. The mass boundary layer thickness $δ_k ≈ d_0 Pe^{1/2}$ is resolved up to 4 cells for each case by increasing the resolution.
**Figure 8(a)** shows that the time evolution of the bubble radius is grid-converged between L9 and L10 for *Sc* = 1. Next, **figure 8(b)** shows that the non-dimensional transfer coefficient is almost constant with time, and *Sh*(t) (non-dimensionalized *k*<sub>L</sub>) reaches a stationary value, while the bubble size is shrinking. The results at *Re*<sub>λ</sub> = 55, *Sc* = 10 confirm this observation, as shown in **figures 8(c,d)**, with convergence obtained between L10 and L11. For *Re*<sub>λ</sub> = 150, *Sc* = 10, the time evolution of the bubble radius and *Sh* is converging to a constant value as resolution is increased to L12. The features in **figure 8** are general to all the cases that we simulated for a wide range of *Re* and *Sc*, as summarized in **table 1**.

We also test the sensitivity of the transfer rate results with the Stefan number, or gas
solubility. The Stefan number is St = α/(1 − α) = 1 for α = 0.5, and St ≈ 0.43 for α =
0.3. The effect of solubility α and concentration of the bubble should not affect the mass
transfer coefficient kL for fully unsaturated cL = 0 conditions as it just rescales the constant
concentration at the bubble boundary kL = (dc/dt)/(αcG − cL). This can also be seen in
**figure 8(b)**, where there is negligible change in kL when α is reduced from 0.5 to 0.3. As
the Reynolds and Schmidt numbers are increased, the boundary layer thickness decreases
as δk ~ Pe−1/2.

The observation of a constant mass transfer coefficient while the bubble size is shrinking contradicts the proposed scaling from Levich (1962) for the mass transfer coefficient in turbulence, *Sh* ∝ (d/L*) *Re*³/₄ *Sc*¹/², which involves a bubble size dependence *d*, as well as our more recent scaling *Sh*/*Sc*¹/² ∝ *Re*¹/² (Farsoiya *et al.* 2021). As the bubble size was constant in Farsoiya *et al.* (2021), the dependency with bubble size had not been tested. We performed additional simulations with the dilute transfer framework from Farsoiya *et al.* (2021), varying the bubble sizes in the same turbulent flow and confirm that the mass transfer coefficient *k*<sub>L</sub> in turbulence does not depend on the bubble size, as demonstrated in **figure 8**. This result can be rationalized by the fact that the smallest length scale of the turbulent flow will control the gas exchange, as described by Lamont & Scott (1970) and the high *Re* regime from Theofanous *et al.* (1976). As will be shown in the next subsection, all the data from the present study and Farsoiya *et al.* (2021) can be described by a high *Re* scaling reading, *Sh*/*Sc*¹/² ∝ *Re*³/⁴, with the Reynolds number and the Sherwood number defined using the integral length scale *L**: *Sh* = *k*<sub>L</sub>*L**/*D*<sub>l</sub> and *Re* = *u*<sub>rms</sub>*L**/*v*<sub>l</sub>, with *u*<sub>rms</sub> the velocity fluctuations.

**Figure 9(a)** shows the concentration of gas transferred in the liquid as a function of time for *Re* = 38, *Sc* = 1. The mass transfer rate dN/dt ∝ −*k*<sub>L</sub>*A* reduces as the bubble size decreases, and reaches zero when the bubble gets dissolved completely as *A* → 0 (figure 9a). The corresponding rate of change of the bubble volume is then dV/dt ∝ dN/dt ∝ −*k*<sub>L</sub>*A*. For a spherical bubble, *V* ∝ *R*(t)³, *A* ∝ *R*(t)², so we have d(*R*(t))³/dt ∝ −*k*<sub>L</sub>*R*(t)², and this implies d*R*/dt ∝ −*k*<sub>L</sub> (or *Sh*). **Figure 9(b)** shows the normalized rate of change of the bubble radius as a function of the Sherwood number, and we indeed observe (d*R*/dt)*L**/(*α*D*<sub>l</sub>) ∝ *Sh*. (The data in **figure 9(b)** have Pearson correlation coefficient ρ((d*R*/dt)*L**/(*α*D*<sub>l</sub>), *Sh*) = 0.94.) As discussed above, the transfer rate *k*<sub>L</sub> and its non-dimensional version *Sh* are constant while the bubble dissolves (see **figures 8b,d**), so the rate of change of the bubble radius is approximately linear with time (which is visible in **figures 8a,c**). Even if the mass transfer coefficient *k*<sub>L</sub> is independent of the bubble size, a larger bubble has a larger interfacial area, leading to a higher mass transfer than a smaller one. The interfacial area is the determining factor for a given *k*<sub>L</sub>, and if a larger bubble breaks into smaller bubbles due to the turbulence, then more surface area is created and the gas transfer increases.
## Bubble-mediated gas transfer

<figure><img src="image_25.png" /><figcaption>Figure 9. (a) Normalized concentration of gas in surrounding liquid (where $c(t_d)$ is the gas concentration in the liquid when the bubble is dissolved completely) for case $Sc = 1$ (highest resolution, L10) for increasing Reynolds number. The gas transfer rate $dc/dt$ decreases and reaches zero as the bubble dissolves fully. (b) Normalized rate of change of the bubble radius $(dR/dt)L^*/(\alpha \mathcal{D}_l)$ as a function of $Sh$ (highest resolution points). A linear scaling is observed.</figcaption></figure>

### 4.4. Model for the mass transfer coefficient

The mass transfer coefficient $k_L$ from the classic penetration theory (Higbie 1935) reads

$$k_L \propto \sqrt{\mathcal{D}_l / \theta}. \qquad (4.2)$$

We observe that the mass transfer coefficient is a function of the bulk turbulence characteristics, but not the bubble size, so we propose the time factor $\theta$ for the gas diffusion in an HIT flow to scale with the Kolmogorov time scale $\tau_\eta = (v_l/\epsilon)^{1/2}$ or integral time scale $\tau_L^* = L^*/u_{rms}$. A similar discussion for scalings of the exchange velocities can be found in Lamont & Scott (1970), Theofanous et al. (1976) and Herlina & Wissink (2016).

Considering the Kolmogorov time scale, the mass transfer coefficient scales as

$$k_L/(v_l\epsilon)^{1/4} \propto Sc^{-1/2}, \qquad (4.3)$$

while considering the integral time scale, the coefficient scales as

$$k_L/u_{rms} \propto Sc^{-1/2} Re^{-1/2}. \qquad (4.4)$$

If we scale $k_L$ by $\mathcal{D}_l/L^*$, and define a turbulent Sherwood number $Sh = k_L L^* / \mathcal{D}_l$, then we obtain

$$Sh/Sc^{1/2} \propto Re^{3/4} \qquad (4.5)$$

and

$$Sh/Sc^{1/2} \propto Re^{1/2}, \qquad (4.6)$$

respectively. The scaling equation (4.5) would be expected to apply only at a high enough Reynolds number while (4.6) would apply to a lower $Re$, similarly to the discussion in Theofanous et al. (1976).
<figure><img src="image_26.png" /><figcaption>Figure 10. (a) Mean non-dimensional mass transfer coefficient Sh after reaching stationary state and R(t) > 0.3R₀, normalized by Sc¹/². The highest (converged) resolution points are shown from table 1. Lines provide the power law for Reⁿ, where Re = L*umorms/vl, with a solid line for n = 3/4, and a dashed line for n = 1/2. The high Re data align with the n = 3/4 exponent (when grid convergence is achieved). Circles, squares and diamonds are dissolving data, while triangles are simulations performed under the dilute approximation; see Farsoiya et al. (2021). All data converge to Sh/Sc¹/² ∝ Re³/⁴. (b) The same scaling but with data plotted as kL/(εv)¹/⁴ as a function of Sc, leading to kL/(εv)¹/⁴ ∝ Sc⁻¹/² (solid line).</figcaption></figure>

Figure 10 shows the turbulent Sherwood number Sh for all the simulations for the cases given in table 1 at the highest numerical resolution (convergence of Sh when resolution is increased is discussed in Appendix A). The error bars depict the minimum and maximum values of Sh after reaching stationary state and before the end of the simulation. All cases with Sc = 1 for all Reynolds numbers (from table 1) are converged to the scaling Sh/Sc¹/² ∝ Re³/⁴. The Sc = 10 cases are converged to the scaling for Re ≤ 400, as the resolution is increased from L10 to L11, while at the highest Sc, the L12 resolution cases are converging. The accepted resolution criteria for direct numerical simulations in the literature (Overholt & Pope 1996; Pope 2001; Schumacher, Sreenivasan & Yeung 2005; Dodd et al. 2021; Farsoiya et al. 2021) are typically kmaxη > 1.5 and kmaxηB > 1.5, where kmax, η and ηB are the maximum resolved wavenumber kmax = πN/L, and the Kolmogorov and Batchelor scales, respectively. The Kolmogorov length scale η = (v³/l/ε)¹/⁴ defines the length scale at which viscous dissipation becomes dominant, while the Batchelor scale is defined as ηB = η/√Sc. In all of the cases, the Kolmogorov length scale is well resolved, with kmaxη > 8, in agreement with the fact that convergence is already achieved in lower resolutions, as shown in figure 5(a). For the highest Péclet number Pe = 4.8 × 10⁴ (Re = 1600, Sc = 30), the Batchelor length scale is resolved up to kmaxηB ≈ 4 (refinement d₀/Δx = 546). The boundary layer thickness is δv ≳ 2η and δk ≳ 2ηB.

The data from dilute gas transfer where bubble size is constant from Farsoiya et al. (2021) also show agreement with (4.5) within narrow error bounds. A new set of simulations for a constant size bubble (Farsoiya et al. 2021) with different diameters
## Bubble-mediated gas transfer

$(d_0/\eta \approx 33, 17, 8, 4)$ is shown by upward triangles at $Re = 800$ in figure 10. This also corroborates that the mass transfer coefficient $k_L$ (or $Sh$) is independent of the bubble size. We note that the constant size dilute mass transfer study used a different numerical method (Farsoiya et al. 2021) and typically higher resolution. At the same resolutions, the error bars between the two methods are comparable. In summary, at the highest resolution and when the Bachelor length scale is resolved, we observe that the transfer coefficient follows the high $Re$ scaling, (4.5) and $Sh/Sc^{1/2} \propto Re^{3/4}$, so that the smallest scales of the turbulent flow drive the mass transfer from the bubble.

The mass transfer scaling can be expressed in terms of $k_L/(\epsilon v_l)^{1/4} \propto Sc^{-1/2}$ and is shown in figure 10(b). A similar conclusion can be made, with convergence to (4.3) at a high enough grid resolution. We note that a higher $Sc$ is achieved for the dilute case and confirms the above discussion.

## 5. Conclusion

We developed a numerical framework to simulate gas transfer across fluid-fluid interfaces including volume change effects. The numerical framework is tested against theoretical and numerical solutions of moving boundary problems. For a bubble rising in a quiescent liquid, the mass transfer coefficient $k_L$ can be described by $k_L = (2/\sqrt{\pi})\sqrt{D_l} U(t)/d(t)$, with $d(t)$ and $U(t)$ the time-varying bubble size and rise velocity. For a bubble in a homogeneous and isotropic turbulence flow, with size initially within the inertial range, we found that the mass transfer coefficient for varying bubble size and Reynolds number is independent of the bubble size. The results are obtained using two types of numerical methods: the one from the present study, which allows for dissolution, and the one presented by Farsoiya et al. (2021) under the dilute approximation, provided that the turbulent microscales and boundary layers are resolved at least up to 4 grid points. The transfer rate in the turbulent flow is controlled by the smallest scale in the flow and can then be described as a function of the Reynolds number and Schmidt number by the scaling $Sh/Sc^{1/2} \propto Re^{3/4}$, where $Sh = k_L L^* / D_l$ is a turbulence Sherwood number. This scaling is equivalent to $k_L \propto Sc^{-1/2}(\epsilon v_l)^{1/4}$, in agreement with the model proposed by Lamont & Scott (1970) and corresponding to the high $Re$ regime from Theofanous et al. (1976). This improved understanding of bubble-mediated gas transfer could be applied to complex systems such as gas transfer related to breaking waves at the air-sea interface or gas transfer in bubble swarms.

**Funding.** This work was supported by the NSF award 2122042 to L.D., the Catalysis Initiative at Princeton, and the Cooperative Institute for Modeling the Earth System (CIMES) between Princeton and NOAA-GFDL. We would like to acknowledge high-performance computing support from the Extreme Science and Engineering Discovery Environment (XSEDE) (Towns et al. 2014), which is supported by the NSF using allocation TG-OCE180010 to L.D. and Cheyenne (doi:10.5065/D6RX99HX) provided by NCAR's Computational and Information Systems Laboratory sponsored by the NSF.

**Declaration of interests.** The authors report no conflict of interest.

**Author ORCIDs.**

[i] Palas Kumar Farsoiya https://orcid.org/0000-0003-4784-3470;

[ii] Quentin Magdelaine https://orcid.org/0000-0002-5705-8313;

[iii] Arnaud Antkowiak https://orcid.org/0000-0002-1534-4290;

[iv] Stéphane Popinet https://orcid.org/0000-0001-9947-297X;

[v] Luc Deike https://orcid.org/0000-0002-4644-9909.
<figure><img src="image_27.png" /><figcaption>Figure 11. Difference between theoretical scaling and numerical data as a function of the grid resolution (level). This is similar to an error convergence plot. We consider $\chi = |Sh - 0.65 Re^{3/4} Sc^{1/2}|/(0.65 Re^{3/4} Sc^{1/2})$ for some cases in **table 1**. The error decreases with increased resolution, converging towards the theoretical scaling. (a) Number of grid cells per initial diameter of the bubble. (b) Number of grid cells per Batchelor scale.</figcaption></figure>

<table><caption>Table 2. Numerical values for Sh/Sc<sup>1/2</sup>.</caption><thead><tr><th rowspan="2">Re<sub>λ</sub></th><th rowspan="2">Sc</th><th colspan="4">Sh/Sc<sup>1/2</sup></th></tr><tr><th>L9</th><th>L10</th><th>L11</th><th>L12</th></tr></thead><tbody><tr><td rowspan="3">38</td><td>1</td><td>21 ± 5</td><td>21 ± 5</td><td>—</td><td>—</td></tr><tr><td>10</td><td>19 ± 10</td><td>20 ± 20</td><td>18 ± 6</td><td>19 ± 0.3</td></tr><tr><td>30</td><td>—</td><td>21 ± 70</td><td>19 ± 6</td><td>22 ± 1</td></tr><tr><td rowspan="3">55</td><td>1</td><td>41 ± 5</td><td>40 ± 5</td><td>—</td><td>—</td></tr><tr><td>10</td><td>47 ± 30</td><td>35 ± 20</td><td>31 ± 15</td><td>42 ± 10</td></tr><tr><td>30</td><td>—</td><td>45 ± 120</td><td>45 ± 30</td><td>48 ± 17</td></tr><tr><td rowspan="3">77</td><td>1</td><td>57 ± 5</td><td>58 ± 5</td><td>—</td><td>—</td></tr><tr><td>10</td><td>94 ± 60</td><td>62 ± 20</td><td>61 ± 10</td><td>57 ± 4</td></tr><tr><td>30</td><td>—</td><td>81 ± 30</td><td>63 ± 7</td><td>58 ± 4</td></tr><tr><td rowspan="3">100</td><td>1</td><td>121 ± 7</td><td>118 ± 2</td><td>—</td><td>—</td></tr><tr><td>10</td><td>163 ± 40</td><td>140 ± 30</td><td>130 ± 15</td><td>118 ± 2</td></tr><tr><td>30</td><td>—</td><td>162 ± 100</td><td>137 ± 30</td><td>133 ± 7</td></tr><tr><td rowspan="3">150</td><td>1</td><td>180 ± 10</td><td>173 ± 7</td><td>—</td><td>—</td></tr><tr><td>10</td><td>496 ± 350</td><td>282 ± 70</td><td>235 ± 70</td><td>186 ± 6</td></tr><tr><td>30</td><td>—</td><td>162 ± 100</td><td>137 ± 30</td><td>238 ± 10</td></tr></tbody></table>

## Appendix A. Grid convergence of the mass transfer rate

Figure 11 shows the normalized difference between the numerical mass transfer coefficient and the theoretical scaling proposed (see (4.5)):

$$ \chi = |Sh - 0.65 Re^{3/4} Sc^{1/2}|/(0.65 Re^{3/4} Sc^{1/2}) \quad (A1) $$

for increasing grid resolutions. The values are $\chi \le 0.5$ and decreasing when going to L11. For higher Péclet numbers, more refined grids are used and allow us to bring the values of $\chi$ below 0.5. Note that the parameter $\chi$ used here is for grid convergence, calculated based on a proportionality constant in (4.5) for best fit in figure 10, and is within $O(1)$ for the data plotted on a log-log scale. Table 2 shows the values for the data points presented in figure 10.
## Bubble-mediated gas transfer

## Appendix B. Computational cost

The computational cost for the present study is given in table 3.

<table><caption>Table 3. Computational cost (with adaptive mesh refinement) of the simulations of bubble dissolution in turbulence.</caption><thead><tr><th>Resolution (2<sup>L</sup>)</th><th>Re<sub>λ</sub></th><th>Re</th><th>Sc</th><th>CPUh</th></tr></thead><tbody><tr><td>L9 ≡ 2<sup>9</sup> to L12 ≡ 2<sup>12</sup></td><td>38</td><td>85</td><td>(1, 10, 30)</td><td>2.04 × 10<sup>5</sup></td></tr><tr><td>L9 ≡ 2<sup>9</sup> to L12 ≡ 2<sup>12</sup></td><td>55</td><td>170</td><td>(1, 10, 30)</td><td>3.43 × 10<sup>5</sup></td></tr><tr><td>L9 ≡ 2<sup>9</sup> to L12 ≡ 2<sup>12</sup></td><td>77</td><td>400</td><td>(1, 10, 30)</td><td>3.27 × 10<sup>5</sup></td></tr><tr><td>L9 ≡ 2<sup>9</sup> to L12 ≡ 2<sup>12</sup></td><td>100</td><td>800</td><td>(1, 10, 30)</td><td>2.11 × 10<sup>5</sup></td></tr><tr><td>L9 ≡ 2<sup>9</sup> to L12 ≡ 2<sup>12</sup></td><td>150</td><td>1600</td><td>(1, 10, 30)</td><td>3.70 × 10<sup>5</sup></td></tr><tr><td>—</td><td>—</td><td>—</td><td>—</td><td>1.45 × 10<sup>6</sup></td></tr></tbody></table>

## REFERENCES

* ALEXIADES, V. & SOLOMON, A.D. 2018 *Mathematical Modeling of Melting and Freezing Processes*. Routledge.

* BERNY, A., DEIKE, L., SÉON, T. & POPINET, S. 2020 Role of all jet drops in mass transfer from bursting bubbles. *Phys. Rev. Fluids* **5** (3), 033605.

* BOTHE, D. & FLECKENSTEIN, S. 2013 A volume-of-fluid-based method for mass transfer processes at fluid particles. *Chem. Engng Sci.* **101**, 283–302.

* BOTHE, D., KOEBE, M., WIELAGE, K., PRÜSS, J. & WARNECKE, H.-J. 2004 Direct numerical simulation of mass transfer between rising gas bubbles and water. In *Bubbly Flows* (ed. M. Sommerfeld), pp. 159–174. Springer.

* BOWMAN, C.W. & JOHNSON, A.I. 1962 Mass transfer from carbon dioxide bubbles rising in water. *Can. J. Chem. Engng* **40** (4), 139–147.

* BRACKBILL, J.U., KOTHE, D.B. & ZEMACH, C. 1992 A continuum method for modeling surface tension. *J. Comput. Phys.* **100** (2), 335–354.

* CANO-LOZANO, J.C., MARTINEZ-BAZAN, C., MAGNAUDET, J. & TCHOUFAG, J. 2016 Paths and wakes of deformable nearly spheroidal rising bubbles close to the transition to path instability. *Phys. Rev. Fluids* **1** (5), 053604.

* CHU, S. & PROSPERETTI, A. 2016 Dissolution and growth of a multicomponent drop in an immiscible liquid. *J. Fluid Mech.* **798**, 787–811.

* CLAASSEN, C.M.Y., ISLAM, S., PETERS, E.A.J.F., DEEN, N.G., KUIPERS, J.A.M. & BALTUSSEN, M.W. 2020 An improved subgrid scale model for front-tracking based simulations of mass transfer from bubbles. *AIChE J.* **66** (4), e16889.

* CLIFT, R., GRACE, J.R. & WEBER, M.E. 2005 *Bubbles, Drops, and Particles*. Courier Corporation.

* COLOMBET, D., LEGENDRE, D., COCKX, A., GUIRAUD, P., RISSO, F., DANIEL, C. & GALINAT, S. 2011 Experimental study of mass transfer in a dense bubble swarm. *Chem. Engng Sci.* **66** (14), 3432–3440.

* COLOMBET, D., LEGENDRE, D., RISSO, F., COCKX, A. & GUIRAUD, P. 2015 Dynamics and mass transfer of rising bubbles in a homogeneous swarm at large gas volume fraction. *J. Fluid Mech.* **763**, 254–285.

* CREPEAU, J. 2007 Josef Stefan: his life and legacy in the thermal sciences. *Expl Therm. Fluid Sci.* **31** (7), 795–803.

* DARMANA, D., DEEN, N.G. & KUIPERS, J.A.M. 2006 Detailed 3D modeling of mass transfer processes in two-phase flows with dynamic interfaces. *Chem. Engng Technol.* **29** (9), 1027–1033.

* DAVIDSON, M.R. & RUDMAN, M. 2002 Volume-of-fluid calculation of heat or mass transfer across deforming interfaces in two-fluid flow. *Numer. Heat Transfer B* **41** (3–4), 291–308.

* DEIKE, L. 2022 Mass transfer at the ocean–atmosphere interface: the role of wave breaking, droplets, and bubbles. *Annu. Rev. Fluid Mech.* **54**, 191–224.

* DEIKE, L. & MELVILLE, W.K. 2018 Gas transfer by breaking waves. *Geophys. Res. Lett.* **45** (19), 10482–10492.
DEISING, D., MARSCHALL, H. & BOTHE, D. 2016 A unified single-field model framework for volume-of-fluid simulations of interfacial species transfer applied to bubbly flows. *Chem. Engng Sci.* **139**, 173–195.

DODD, M. 2017 Direct numerical simulation of droplet-laden isotropic turbulence. PhD thesis, University of Washington, Seattle, WA.

DODD, M.S. & FERRANTE, A. 2016 On the interaction of Taylor length scale size droplets and isotropic turbulence. *J. Fluid Mech.* **806**, 356.

DODD, M.S., MOHADDES, D., FERRANTE, A. & IHME, M. 2021 Analysis of droplet evaporation in isotropic turbulence through droplet-resolved DNS. *Intl J. Heat Mass Transfer* **172**, 121157.

DUNCAN, P.B. & NEEDHAM, D. 2004 Test of the Epstein-Plesset model for gas microparticle dissolution in aqueous media: effect of surface tension and gas undersaturation in solution. *Langmuir* **20** (7), 2567–2578.

ELGHOBASHI, S. 2019 Direct numerical simulation of turbulent flows laden with droplets or bubbles. *Annu. Rev. Fluid Mech.* **51** (1), 217–244.

EPSTEIN, P.S. & PLESSET, M.S. 1950 On the stability of gas bubbles in liquid-gas solutions. *J. Chem. Phys.* **18** (11), 1505–1509.

FARSOIYA, P., POPINET, S. & DEIKE, L. 2022a Stefan problem. http://basilisk.fr/sandbox/farsoiya/phase_change/stefan Problem.c.

FARSOIYA, P., POPINET, S. & DEIKE, L. 2022b Epstein-Plesset test. http://basilisk.fr/sandbox/farsoiya/phase_change/epstein-plesset.c.

FARSOIYA, P.K., POPINET, S. & DEIKE, L. 2021 Bubble-mediated transfer of dilute gas in turbulence. *J. Fluid Mech.* **920**, A34.

FLECKENSTEIN, S. & BOTHE, D. 2015 A volume-of-fluid-based numerical method for multi-component mass transfer with local volume changes. *J. Comput. Phys.* **301**, 35–58.

FORTESCUE, G.E. & PEARSON, J.R.A. 1967 On gas absorption into a turbulent liquid. *Chem. Engng Sci.* **22** (9), 1163–1176.

GARBE, C.S., et al. 2014 Transfer across the air-sea interface. In *Ocean-Atmosphere Interactions of Gases and Particles* (ed. P.S. Liss & M.T. Johnson), pp. 55–112. Springer.

GENNARI, G., JEFFERSON-LOVEDAY, R. & PICKERING, S.J. 2022 A phase-change model for diffusion-driven mass transfer problems in incompressible two-phase flows. *Chem. Engng Sci.* **259**, 117791.

GOTOH, T., YEUNG, P.K., DAVIDSON, P.A., KANEDA, Y. & SREENIVASAN, K.R. 2013 Passive scalar transport in turbulence: a computational perspective. In *Ten Chapters in Turbulence* (ed. P.A. Davidson, Y. Kaneda & K.R. Sreenivasan), pp. 87–131. Cambridge University Press.

GUMULYA, M., UTIKAR, R.P., PAREEK, V.K., EVANS, G.M. & JOSHI, J.B. 2020 Dynamics of bubbles rising in pseudo-2D bubble column: effect of confinement and inertia. *Chem. Engng J.* **405**, 126615.

HAROUN, Y., LEGENDRE, D. & RAYNAL, L. 2010 Volume of fluid method for interfacial reactive mass transfer: application to stable liquid film. *Chem. Engng Sci.* **65** (10), 2896–2909.

HERLINA, H. & WISSINK, J.G. 2016 Isotropic-turbulence-induced mass transfer across a severely contaminated water surface. *J. Fluid Mech.* **797**, 665–682.

HERLINA, H. & WISSINK, J.G. 2019 Simulation of air-water interfacial mass transfer driven by high-intensity isotropic turbulence. *J. Fluid Mech.* **860**, 419–440.

HIGBIE, R. 1935 The rate of absorption of a pure gas into a still liquid during short periods of exposure. *Trans. AIChE* **31**, 365–389.

VAN HOOFT, J.A., POPINET, S., VAN HEERWAARDEN, C.C., VAN DER LINDEN, S.J.A., DE ROODE, S.R. & VAN DE WIEL, B.J.H. 2018 Towards adaptive grids for atmospheric boundary-layer simulations. *Boundary-Layer Meteorol.* **167** (3), 421–443.

KEELING, R.F. 1993 On the role of large bubbles in air-sea gas exchange and supersaturation in the ocean. *J. Mar. Res.* **51** (2), 237–271.

LAMONT, J.C. & SCOTT, D.S. 1970 An eddy cell model of mass transfer into the surface of a turbulent liquid. *AIChE J.* **16** (4), 513–519.

LEGENDRE, D. & MAGNAUDET, J. 1999 Effet de l'accélération d'un écoulement sur le transfert de masse ou de chaleur à la surface d'une bulle sphérique. *C. R. l'Acadé. Sci.* **327** (1), 63–70.

LEVICH, V.G. 1962 *Physicochemical hydrodynamics*. Prentice-Hall.

LIANG, J.-H., MCWILLIAMS, J.C., SULLIVAN, P.P. & BASCHEK, B. 2011 Modeling bubbles and dissolved gases in the ocean. *J. Geophys. Res.* **116**, C03015.

LOHSE, D. & ZHANG, X., et al. 2015 Pinning and gas oversaturation imply stable single surface nanobubbles. *Phys. Rev. E* **91** (3), 031003.

LOISY, A. & NASO, A. 2017 Interaction between a large buoyant bubble and turbulence. *Phys. Rev. Fluids* **2** (1), 014606.
LOISY, A., NASO, A. & SPELT, P.D.M. 2017 Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction. *J. Fluid Mech.* **816**, 94–141.

LÓPEZ-HERRERA, J.M., GANAN-CALVO, A.M., POPINET, S. & HERRADA, M.A. 2015 Electrokinetic effects in the breakup of electrified jets: a volume-of-fluid numerical study. *Intl J. Multiphase Flow* **71**, 14–22.

LUNDGREN, T.S. 2003 Linearly forced isotropic turbulence. In *Annual Research Briefs CTR* (ed. P. Moin & N.N. Mansour), pp. 461–473. Stanford University.

MAES, J. & SOULAINE, C. 2020 A unified single-field volume-of-fluid-based formulation for multi-component interfacial transfer with local volume changes. *J. Comput. Phys.* **402**, 109024.

MAGDELAINE, Q. 2019 Hydrodynamique des films liquides hétérogènes. PhD Thesis, Sorbonne Université. https://tel.archives-ouvertes.fr/tel-03139815.

MALAN, L. 2018 Direct numerical simulation of free-surface and interfacial flow using the VOF method: cavitating bubble clouds and phase change. PhD thesis, University of Cape Town.

MARIĆ, T., KOTHE, D.B. & BOTHE, D. 2020 Unstructured un-split geometrical volume-of-fluid methods – a review. *J. Comput. Phys.* **420**, 109695.

MARSCHALL, H., HINTERBERGER, K., SCHÜLER, C., HABLA, F. & HINRICHSEN, O. 2012 Numerical simulation of species transfer across fluid interfaces in free-surface flows using OpenFOAM. *Chem. Engng Sci.* **78**, 111–127.

MAXWORTHY, T., GNANN, C., KÜRTEN, M. & DURST, F. 1996 Experiments on the rise of air bubbles in clean viscous liquids. *J. Fluid Mech.* **321**, 421–441.

MOORE, D.W. 1965 The velocity of rise of distorted gas bubbles in a liquid of small viscosity. *J. Fluid Mech.* **23** (4), 749–766.

MOSTERT, W., POPINET, S. & DEIKE, L. 2021 High-resolution direct simulation of deep water breaking waves: transition to turbulence, bubbles and droplet production. arXiv:2103.05851.

NASO, A. & PROSPERETTI, A. 2010 The interaction between a solid particle and a turbulent flow. *New J. Phys.* **12** (3), 033040.

OVERHOLT, M.R. & POPE, S.B. 1996 Direct numerical simulation of a passive scalar with imposed mean gradient in isotropic turbulence. *Phys. Fluids* **8** (11), 3128–3148.

PEÑAS-LÓPEZ, P., PARRALES, M.A., RODRÍGUEZ-RODRÍGUEZ, J. & VAN DER MEER, D. 2016 The history effect in bubble growth and dissolution. Part 1. Theory. *J. Fluid Mech.* **800**, 180–212.

PERRARD, S., RIVIÈRE, A., MOSTERT, W. & DEIKE, L. 2021 Bubble deformation by a turbulent flow. *J. Fluid Mech.* **920**, A15.

PESKIN, C.S. 1972 Flow patterns around heart valves: a numerical method. *J. Comput. Phys.* **10** (2), 252–271.

PINELLI, M., HERLINA, H., WISSINK, J.G. & UHLMANN, M. 2022 Direct numerical simulation of turbulent mass transfer at the surface of an open channel flow. *J. Fluid Mech.* **933**, A49.

POPE, S.B. 2001 *Turbulent flows*. Cambridge University Press.

POPINET, S. 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. *J. Comput. Phys.* **228** (16), 5838–5866.

POPINET, S. 2015 A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations. *J. Comput. Phys.* **302**, 336–358.

POPINET, S. 2018 Numerical models of surface tension. *Annu. Rev. Fluid Mech.* **50**, 49–75.

POPINET, S. & COLLABORATORS 2013–2022 Basilisk. http://basilisk.fr.

REICHL, B. & DEIKE, L. 2020 Contribution of sea-state dependent bubbles to air-sea carbon dioxide fluxes. *Geophys. Res. Lett.* **47**, e2020GL087267.

RISSO, F. 2018 Agitation, mixing, and transfers induced by bubbles. *Annu. Rev. Fluid Mech.* **50**, 25–48.

RIVIÈRE, A., MOSTERT, W., PERRARD, S. & DEIKE, L. 2021 Sub-Hinze scale bubble production in turbulent bubble break-up. *J. Fluid Mech.* **917**, A40.

ROGHAIR, I. 2012 Direct numerical simulations of hydrodynamics and mass transfer in dense bubbly flows. PhD thesis, Eindhoven University of Technology.

ROSALES, C. & MENEVEAU, C. 2005 Linear forcing in numerical simulations of isotropic turbulence: physical space implementations and convergence properties. *Phys. Fluids* **17** (9), 095106.

SATO, T., JUNG, R.-T. & ABE, S. 2000 Direct simulation of droplet flow with mass transfer at interface. *J. Fluids Engng* **122** (3), 510–516.

SCAPIN, N., COSTA, P. & BRANDT, L. 2020 A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows. *J. Comput. Phys.* **407**, 109251.

SCARDOVELLI, R. & ZALESKI, S. 1999 Direct numerical simulation of free-surface and interfacial flow. *Annu. Rev. Fluid Mech.* **31** (1), 567–603.

SCHUMACHER, J., SREENIVASAN, K.R. & YEUNG, P.K. 2005 Very fine structures in scalar mixing. *J. Fluid Mech.* **531**, 113–122.
SIDEMAN, S. 1966 The equivalence of the penetration and potential flow theories. *Ind. Engng Chem.* **58** (2), 54–58.

STEFAN, J. 1891 Über die theorie der eisbildung, insbesondere über die eisbildung im polarmeere. *Ann. Phys.* **278** (2), 269–286.

TAKEMURA, F. & YABE, A. 1998 Gas dissolution process of spherical rising gas bubbles. *Chem. Engng Sci.* **53** (15), 2691–2699.

TANGUY, S., SAGAN, M., LALANNE, B., COUDERC, F. & COLIN, C. 2014 Benchmarks and numerical methods for the simulation of boiling flows. *J. Comput. Phys.* **264**, 1–22.

THEOFANOUS, T.G., HOUZE, R.N. & BRUMFIELD, L.K. 1976 Turbulent mass transfer at free, gas-liquid interfaces, with applications to open-channel, bubble and jet flows. *Intl J. Heat Mass Transfer* **19** (6), 613–624.

TOWNS, J., et al. 2014 Xsede: accelerating scientific discovery. *Comput. Sci. Engng* **16** (5), 62–74.

TRYGGVASON, G., SCARDOVELLI, R. & ZALESKI, S. 2011 *Direct Numerical Simulations of Gas–Liquid Multiphase Flows*. Cambridge University Press.

WEINBERG, M.C. & SUBRAMANIAN, R.S. 1980 Dissolution of multicomponent bubbles. *J. Am. Ceram. Soc.* **63** (9–10), 527–531.

WEINER, A. & BOTHE, D. 2017 Advanced subgrid-scale modeling for convection-dominated species transport at fluid interfaces with application to mass transfer from rising bubbles. *J. Comput. Phys.* **347**, 261–289.

WOOLF, D.K. & THORPE, S.A. 1991 Bubbles and the air–sea exchange of gases in near-saturation conditions. *J. Mar. Res.* **49** (3), 435–466.

ZHU, X., VERZICCO, R., ZHANG, X. & LOHSE, D. 2018 Diffusive interaction of multiple surface nanobubbles: shrinkage, growth, and coarsening. *Soft Matt.* **14** (11), 2006–2014.