Experimental and numerical investigation of reactive species transport around a small rising bubble

In this article, we present experimental and numerical techniques to investigate the transfer, transport, and reaction of a chemical species in the vicinity of rising bubbles. In the experiment, single oxygen bubbles of diameter d b = 0.55 . . . 0.85 mm are released into a measurement cell filled with tap water. The oxygen dissolves and reacts with sulfite to sulfate. Laser-induced fluorescence is used to visualize the oxygen concentration in the bubble wake from which the global mass transfer coefficient can be calculated. The ruthenium-based fluorescent dye seems to be surface active, such that the rise velocity is reduced by up to 50 % compared to the experiment without fluorescent dye and a recirculation zone forms in the bubble wake. To access the local mass transfer at the interface, we perform complementary numerical simulations. Since the fluorescence tracer is essential for the experimental method, the effect of surface contamination is also considered in the simulation. We employ several improvements in the experimental and numerical procedures which allow for a quantitative comparison (locally and globally). Rise velocity and mass transfer coefficient agree within a few percents between experiment, simulation and literature results. Because the fluorescence tracer is frequently used in mass transfer experiments, we discuss its potential surface activity.


Introduction
The mass transfer from a gaseous, dispersed phase into the surrounding liquid phase is a major task within the chemical industry as well as in bio-, food-or environmental engineering. The effective utilization of the gaseous phase with a minimum in energy consumption, optimal contact times and high mass transfer rates has gained more and more importance for process intensification and optimization. Nevertheless, a reliable and exact design of multiphase reactors is one of the unsolved challenges in process engineering since local mass transfer processes in dense bubbly flows cannot be investigated easily. Instead, it is common practice to investigate the behavior of a single bubble and then to apply case-dependent correction factors to account for the effect of bubble-bubble interactions. But also correlations for single fluid particles are mainly empirical or semi-empirical and often ambiguous. Semi-empirical correlations for the Sherwood number (Sh) typically have the form of equation (1) below. The idea behind the equation is as follows: the transport of a passive chemical species around a fluid particle rising in a stagnant liquid may be characterized by the Reynolds (Re = d b U b /ν) and Schmidt (Sc = D/ν) number, where d b is the bubble diameter, U b is the rise velocity, ν is the liquid side kinematic viscosity, and D is the molecular diffusivity of the dissolved gas. A simplified analysis based on boundary layer theory leads to a scaling of Sh ∼ P e 1/2 with the Pélect number defined as P e = ReSc (see for example reference [16]). The actual scaling may deviate, for example, if the flow around the bubble is complex (recirculation, vortex shedding), if the bubble is deformed, or if the shape oscillates (wobbling). The simple scaling law also breaks down if the transport is not convection-dominated. For example, as Re → 0, the Sherwood number of a spherical particle should approach Sh = 2 and not zero. This theoretical considerations presumably led Brauer [4] to suggest the empirically corrected equation containing case-specific constants C 1 , a, and b. Another critical factor is the presence of surface contamination. As contaminant, we consider any substance in the system which adsorbs at the interface modifying the surface properties, in particular, the surface tension of the gas-liquid interface. Even the presence of tiny amounts of such a substance may significantly modify the mass transfer performance of a system because the interface is partially immobilized. From this point of view, most systems are influenced by contamination. Already the usage of tap water instead of purified water may reduce the rise velocity up to 50 % (see figure 7.3 [5, p. 172]). Thus, it is not too surprising that experimental results often agree well with the scaling behavior of solid particles Sh ∼ Re 1/2 Sc 1/3 . A frequently cited correlation is that of Frössling [7]: Note that the coefficient preceding the Schmidt number in equation (2) is actually based on experimental results of heat transfer at falling droplets. If the transfer species reacts close to the interface with some bulk component, the (transfer) species boundary layer is depleted and the mass transfer enhanced. Instead of defining correlations for a reactive Sherwood number, it is more common to apply an enhancement factor to the non-reactive value, e.g., Sh r = Sh E. A frequently used correlation for E is based on film theory and, for slow reactions, reads as where the Hatta number is defined as Ha = √ kD/k L for a first-order reaction with the reaction rate k and the non-reactive mass transfer coefficient k L .
Even though experimental and numerical results reflect theoretical scaling laws, there is some ambiguity when it comes to the precise quantitative prediction of mass transfer, or when results differ from the expected scaling behavior. Different measurement or simulation techniques will lead to different correction factors that match the data. The purity of materials used in the experiment or the cleanliness of apertures may significantly affect the obtained results. For example, it is known that rise path and velocity of bubbles rising in water depend strongly on the water's purity. The purity of standard tap water may change on a daily basis in one location, and it certainly varies between different locations (cities). Many times it is not possible or not even wanted to exclude such effects because the ultimate goal is to predict the performance of real-world processes, and, therefore, experiments should resemble their industrial counterparts as closely as possible. On the other hand, experiments also aim to focus on isolated aspects to gain an understanding of the underlying mechanisms. The gap between these two conflicting goals can be reduced by the aid of numerical tools. Simulations are based on mathematical models which aim to describe physical phenomena. When simulation and experiment agree within a few percents with respect to the main performance parameters like rise velocity or mass transfer coefficient, the dominating effects in the experiment can be considered as understood. Furthermore, simulations contain a higher degree of information which allows to comprehend why predictions agree or disagree with expectations. To leverage this synergy effect between experiment and simulation, numerical tools have to include all the principal physical phenomena and they have to be able to investigate these effects on all relevant time and length scales.
The present work focuses on a cross-validation between experiments based on the (planar) laser-induced fluorescence (LIF) technique [3,6,13,29,14,11,9] and numerical results based on Arbitrary Eulerian-Lagrangian (ALE) Interface-Tracking [10,26,27,18,19,20]. Both techniques allow to record local transport phenomena with high temporal and spatial resolution. We reproduced the experimental setting by Kück et al. [29,14], in which relatively small single rising oxygen bubbles were investigated. The bubble diameter ranges between 0.55 mm and 0.85 mm, such that a mostly rectilinear path is attained. The aim is to reduce the complexity by avoiding the additional influence of path instability and shape oscillations. The reaction system is the same as in our previous works, too. The sulfite-sulfate oxidation has several advantages and disadvantages, which are briefly discussed in section 2. Thanks to developments regarding the frame rate and the laser sheet structure on the experimental side, and a numerical technique that is able to handle surface contamination as well as realistic material properties, we were able to find a promising qualitative and quantitative agreement between experiment, simulation and literature results. The investigation techniques and their improvements are explained in section 3. Results for the rise velocity, the mass transfer coefficient and the oxygen distribution in the bubble wake are presented in section 4. The latter section also discusses the ruthenium-based fluorescence tracer as a potential surface-active agent (surfactant) and shows how even a small amount of surfactant drastically changes the mass transfer characteristics.

Sulfite-sulfate-reaction
Understanding the mixing in gas-liquid contactors is a crucial step for process-intensification. A possible strategy to improve the mixing without additional energy consumption is to exploit the bubble wake dynamics as a local "mixing device". To investigate and quantify this strategy, the reaction timescale must be tunable to control precisely where the reaction takes place. The selection of a suitable reaction system is even more challenging for mass transfer investigations based on laser-induced fluorescence techniques, because the transfer species must be fluorescent or quenching. For a gas-liquid interface, a characteristic timescale describing the contact time may be defined as τ conv = d b /U b . The timescale of a first-order reaction is simply the inverse of the reaction rate constant τ r = 1/k. An important dimensionless group defined as the ratio of both timescales is the Damköhler number Da = τ conv /τ r = kd b /U b . For small values of Da, i.e. Da 1, the reaction takes place distributed over significant portions of the liquid bulk. The mass transfer enhancement will be negligible, and the influence of local bubble hydrodynamics on the reaction will be small. For small reaction timescales, i.e. Da 1, the reaction takes place in the immediate vicinity of the interface. The enhancement will be strong, but it won't be influenced much by the flow in the bubble wake. Additional criteria for a suitable reaction system for method validation may be summarized as follows: 1) The physisorption of the transfer species must be known.
2) At least one of the educts or products has to be fluorescent or quenching.
3) The reaction path should be well-defined and known.
4) The reaction timescale should be tunable such that Da ≈ 1. For the present work, we chose the oxidation of sodium sulfite to sulfate, formally written as While the overall reaction is very simple, the detailed reaction mechanism is rather complex due to radical reactions [2]. Actually, there is still some controversy in the literature about the exact reaction steps [21,22,15,12]. Despite the fact that not all details about the reaction kinetics are clarified yet, we chose this system because it has been widely used in previous studies of reactive mass transfer in different groups [13,29,14,24] and because our aim is to show the progress made in the experimental and numerical techniques.

Methods of investigation
3.1. Experimental analysis.
3.1.1. Experimental setup. The accurate investigation of mass transfer at small bubbles rising freely on a rectilinear path requires an experimental setup and procedure that allows the generation of bubbles with reproducible shape and trajectories. The generation of rectilinear rising bubbles with a size smaller than 1 mm was already performed by Ohl [17] and was used previously within the work of Kück et al. [13,29,14]. The technique uses two injection valves, one to control the flow rate of the liquid and gaseous flow, and one to adjust the bubble volume by setting the opening period via a function generator; see figure 1 a). Additionally, the pressure has to be controlled very accurately with precision pressure regulators (range 0.05 to 2 bar, adjustable in 0.001 bar steps). With this technique bubbles from 0.3 to 3 mm can be generated [17]. Figure 1 c) shows the flow scheme of the basic setup. The setup was used in a slightly different configuration in Timmermann et al. [24]. Heart of the setup is the measuring cell (cross section 150 × 150 mm 2 ) which consists of four glass walls with bottom and lid. Both base and top are made of stainless steel to allow different reaction systems. The bubble generator is placed in the center of the bottom part.  For the investigation of local mass transfer at rising gas bubbles, laser-induced fluorescence (LIF) is applied to visualize oxygen concentration fields. The oxygen sensitive dye Dichlorotris-(1,10-phenanthroline)-ruthenium(II)hydrate (c Ru = 30 mg/L; by Sigma-Aldrich ) is used. The fluorescence of this dye shows a dependence on the oxygen concentration and can be described with a Stern-Volmer correlation [29]. The Stern-Volmer constant and the fluorescent intensity in the absence of oxygen are determined by recording grey level images for several different predefined uniform oxygen concentrations. The oxygen is measured with a PreSens oxygen sensor. Additionally, the dye solution contains cobalt(II)sulfate (c Co = 16.67 mg/L) as a catalyst for the sodium sulfite reaction. Excitation of the fluorescent dye is performed with a Nd:YLF laser (wavelength 527 nm, pulse width > 210 ns, pulse repetition rate 1 kHz, Continuum ). The laser beam is widened with light sheet optics to obtain a planar laser sheet in the direction of the rectilinear bubble rise. The emitted fluorescence light is recorded with a PCO Dmiax HS2 (1000 f ps), protected from direct laser radiation by a bandpass filter perpendicular to the laser sheet (ILA 5150 GmbH, center wavelength 590 ± 2 nm, half-power bandwidth 20 ± 2 nm, transmission > 84 %). The camera is equipped with an Infinity K2/Distamax objective to achieve high magnifications (field of view: 5.9 × 7.6 mm). The experimental arrangement for the LIF measurements is shown in figure 2.
To visualize the mass transfer in case of physical and reactive mass transfer, it is necessary to work under oxygen-free conditions. Hence the setup is purged with nitrogen for 30 min, and simultaneously the dye solution is filled into the desorption column and saturated with nitrogen. The oxygen concentration is monitored during the whole measurement within the desorption column and measuring cell. If the oxygen is depleted and the setup is fully purged with nitrogen, the dye solution is filled into the measuring cell.
3.2. Image processing. Due to a non-uniform illumination of the recorded images, a background correction is necessary to obtain reliable concentration field information. The performed correction is based on the work of Dani et al. [6]. First a background image sequence in oxygen desorbed dye solution is recorded. From this sequence, an averaged image is computed and the recorded raw images, see figure 3 a), of the rectilinear rising bubble are divided by the averaged image. Within this corrected image, a nearly uniform background could be obtained; see figure 3 b).
In contrast to the experiments of Kück et al. [29,14], high-speed recordings of the bubble rise were performed to obtain averaged images with small noise information. Therefore the background  corrected images were automatically processed with edge detection in Matlab to detect the bubble center on each image based on the reflection ring inside the bubble. The found bubble center was used to define a region of interest (ROI)-frame of 400 × 800 pixels such that the bubble center is always at a position of 200 pixels in x-and 150 pixels in y-direction. Afterwards, the ROI images were averaged to obtain an image with reduced noise level; see figure 3 c). Nevertheless, sub-pixel shifts of the bubble center could not be compensated with this technique, so that a small blur of the bubble and its wake structure also occurs.

Numerical analysis.
3.3.1. Numerical setup. We used an enhanced version of the bubbleInterTrackFoam [10,26,27] solver contained in version 3.1 of foam-extend 1 (community-driven version of OpenFOAM R ) to perform numerical simulations of clean and contaminated bubbles. The interface-tracking method is complex and computationally expensive, but currently also the most accurate approach to simulate multiphase flows of single fluid particles. One momentum equation is solved for each bulk phase, and the exchange of momentum between the phases occurs via a direct implementation of the jump conditions at the interface. The exact fulfillment of the jump conditions is a major advantage over all methods solving a one-field formulation of the momentum equation; e.g. volume-of-fluid, level-set, front-tracking approaches. Furthermore, the explicit representation of the interface by the computational grid allows to solve surface transport equations. On the other hand, interface tracking solvers are expensive (iterative inter-phase coupling, mesh motion) and less robust than other approaches. The main modifications are a sorption-library [18,19] to simulate the effect of soluble surfactant and a subgrid-scale model to approximate the convection-dominated transport of surfactant in the boundary layer forming at the bubble interface [20]. The sorption-library contains models for fast (diffusion-controlled) and slow (kinetically controlled) sorption as well as corresponding models for the sorption-isotherm. One of the main challenges in this study is the presence of surface-active agents (surfactants) in the experimental system, which alter the properties of the gas-liquid interface. The contaminant source is very likely related to the addition of the fluorescent dye, as will be discussed in section 4.4. From the experimental results of the bubble rise velocity presented in section 4.1, it can be assumed that there is at least one surfactant present in the system. We do not know the exact source of the contamination, however, there is one general observation to our advantage: if the surface contamination is sufficiently high, the bubble reaches a steady velocity which is, over a wide range of surfactant concentrations, independent from the exact concentration value in the liquid bulk [20]. The bubble dynamics in the initial transient stage may vary from case to case, but here we only consider the steady-state regime. Even though the steady-state rise velocity may not depend on the surfactant concentration, it could still depend on surfactant-specific material parameters like the molecular diffusivity of the surfactant in the liquid bulk and on the interface, or the Langmuir constant. Therefore, we tested the behavior under different bulk concentrations of two commonly used surfactants, namely Triton X-100 and C 12 DMPO. In the tested parameter range, the rise velocity never changed more than ±3 % around the mean of all test-cases. This observation is supported by the fact that in most of the literature concerned with the rise of single bubbles in a contaminated liquid (see for example [25,1,20]) the hydrodynamic drag is comparable to the one of a rising solid particle. No strong dependency on the specific surfactant is observed or included in correlations for the drag coefficient. Of course, these assumptions may not hold true in general, but the results presented in the sections hereafter suggest that, at least for the present study, these arguments seem to be justified. For completeness, the exact surfactant properties used in the simulations of contaminated bubbles are reported in table 1. The initial surfactant bulk concentration was set to c 0 = 0.05 mol/m 3 , a very high value which causes the bubble to reach quickly a "fully contaminated" state (a state in which the Marangoni forces remain approximately constant). For a more comprehensive discussion of the properties in table 1 and the sorption model, the reader is referred to the sections 2-4 in [20]. Table 1. Surfactant (C 12 DMPO) properties, fast Langmuir adsorption model parameters. c Σ ∞ -saturated surface concentration, a L -Langmuir constant, D -molecular diffusivity in the liquid bulk, D Σ -molecular diffusivity on the interface, Ttemperature Further simulation parameters are the bulk phase properties of liquid (water) and gas (oxygen) compiled in table 2. Note that in the numerical model the dissolved oxygen is treated as a passive scalar meaning that it does not influence the liquid properties. Moreover, the volume change of the bubble due to the oxygen is neglected. The oxygen concentration fields presented in later chapters result from the solution of a convection-diffusion equation in the liquid phase. A constant value of c| Σ,O 2 = 1.331 mol/m 3 (the saturation concentration of oxygen in water) was set as a boundary condition at the interface. The molecular diffusivity of oxygen D O 2 = 2.0 · 10 −9 m 2 /s is by a factor of four larger than that of the surfactant, however, the concentration boundary layer is still extremely thin such that a subgrid-scale model is necessary to approximate the oxygen transfer accurately. The basic modeling idea is unchanged compared to previous works [31,20]. The analytical solution of a substitute problem is used to improve the numerical solution of the species transport equation close to the interface, where the computational mesh is not fine enough to resolve the steep concentration profile. This approach was used previously for physical species transfer [31] and for the sorption of surfactant [20]. To handle also species transfer with a firstorder chemical reaction, we derived modified model equations based on a "reactive" version of the substitute problem, which can be found in A.3. Table 2. Hydrodynamic material properties for the liquid (+) and gas (−) phase. ρ -density, µ -dynamic viscosity, σ -surface tension (clean).
Due to the enormous computational costs of the numerical method, we simulated only three different bubble diameters d b = 0.70/0.75/0.80 mm. After the bubbles reached a steady-state, velocity field and bubble shape were frozen, and the reactive mass transfer was solved. This procedure allows to investigate a larger range of reaction rates since the solution of a single transport equation requires significantly less computational effort. The meshes have about 220/241/278 thousand bulk and 2850/3124/3607 interface cells for the bubble diameters d b = 0.70/0.75/0.80 mm, respectively. Figure 4 shows, as an example, the mesh for d b = 0.75 mm. We also ran simulations on meshes with approximately twice as many bulk cells and found the rise velocity to vary by less than 0.2 % for both the clean and contaminated case.

3.3.2.
Post-Processing of concentration fields. The measuring principle of fluorescence quenching is based on the effect that oxygen molecules absorb energy of the excited fluorescence tracer. As a result, the fluorescence intensity in regions rich of oxygen is reduced. Based on a calibration curve, gray shades in the recorded camera images can be related to quantitative oxygen concentration values. The laser sheet excites the tracer molecules in a plane parallel to the rise direction and therefore allows to visualize the concentration field in the bubble wake. A similar result may be extracted from the simulation by "slicing" through the computational domain. Mathematically, such a slice or plane is a 2D object with zero thickness. For the comparison of local fields between experiment and simulation, it is essential to understand how far the laser-sheet represents such an object. Figure 5 is a preview on the results presented in section 4.2. The shape of the wake and the angle at which the flow detaches from the interface agree very well between experiment and simulation. However, there are some areas of higher concentration which appear in the experimental recording but not in the simulation. A thorough discussion follows in section 4.2. The laser-sheet has a width of approximately 0.5 mm in the region where it intersects the rising bubble. When considering bubble diameters of less than 1.0 mm this thickness is not negligible. The numerical post-processing step introduced below is a simple attempt to investigate the influence of the finite laser-sheet width. It is conceivable that the measured concentration field is a depth-averaged version of the 3D field in the bubble wake, where the integration/averaging length corresponds to the width of the laser-sheet. Note that this will not alter the absolute concentration value since the calibration curve already contains the same effect. It is more likely that the observed concentration field might become blurry or smeared out. To verify or falsify this assumption, we emulate the same effect in the post-processing of the numerical results. Figure 6 depicts how (triangulated) slices of the concentration field are extracted. This intermediate step is necessary because the computational mesh does not have a fixed underlying structure (like a Cartesian mesh) such that it is not straight-forward to integrate or average the concentration field in a given cuboid. With the extracted slices it is possible to reconstruct the field in the volume in between. A simple linear profile between two x-y-plane at z i−1 and z i leads to the formula where z may be thought of as a coordinate perpendicular to the laser-sheet. The average over N sec sections can then be computed using the trapezoidal rule where z i−1 and z i are positioned in the planes at the start and end of each section, respectively.
The variables x and y in equation (6) were dropped for simplicity. For an equidistant spacing between the planes, the formula reduces to 4. Results 4.1. Rise velocity. As mentioned above, the fluorescence tracer seems to reduce the mobility of the bubble surface. Figure 7 shows the rise velocities obtained from our experiments and simulations compared to reference values from the literature. Relationships for the rise velocity are usually expressed in dimensionless form using the drag coefficient c D . Equating drag and buoyancy force allows to compute a steady-state velocity; see, e.g. equation (5) in [25]): In general, U b cannot be calculated directly from expression (8), unless a simple correlation for c D is used 2 .    the fluorescent dye are close to Tomiyama's reference curve (based on experiments). Note that the results without the fluorescence tracer stem from preliminary tests conducted before the mass transfer experiments. For these cases, rise velocity and diameter are estimated based on edgedetection such that the error margins were small and are therefore not depicted in the plot. On the other hand, in the mass transfer experiments based on LIF, there is only a low gray level difference between the bubble surface and the surrounding liquid. Furthermore, the gray level varies due to the bubble's wake and shadow. Therefore, it is not possible to apply the edgedetection algorithm. Thus, for these cases, the bubble is assumed to be spherical, and velocity and size are determined based on the bubble's shadow; see figure 9. This measurement principle is less accurate and more noisy; see table 3 for a quantification. Table 3. Minimum, maximum and mean relative standard deviation for the measured U b and d b .
After adding the fluorescent dye to the system, the observed terminal velocities are significantly reduced and agree well with the reference for fully contaminated systems. The numerically predicted rise velocities for contaminated bubbles are close to the reference curve for solid particles (or fully contaminated bubbles), too. Figure 8 shows the internal and external fluid flow predicted by the simulations. The interface velocity is significantly reduced and almost stagnant (in a reference frame moving with the bubble center). Furthermore, a recirculation zone forms in the bubble wake. The local Marangoni forces associated with the fully contaminated state are distributed over the entire upper bubble hemisphere. A more detailed discussion of the forces acting on the bubble can be found in [20].

4.2.
Oxygen concentration in the wake. Figure 9 shows the averaged ROI-images in the experiment for several diameters with and without superimposed sodium sulfite reaction. With increasing diameter, the region rich of oxygen becomes broader and longer, as previously reported by Wasowski and Blaß [30]. The observed width of the wake structure indicates the presence of a stagnant ring vortex below the bubble, which is in agreement with the streamlines depicted in figure 8. When sodium sulfite is added to the experiment, the reaction consumes oxygen. Note that the sulfite concentration is at least one order of magnitude higher than that of oxygen, which is why we discuss the addition of sulfite like an increase of the reaction rate in a pseudo-first-order reaction. For the slowest reaction investigated, a small oxygen-drained region in the wake can be observed, which corresponds to the core of the ring-vortex. The transport of additional oxygen into the stagnant recirculation zone happens mainly due to diffusion and is therefore rather slow, such that the oxygen is completely consumed by the reaction. Moreover, the oxygen distribution is more confined compared to the non-reactive case. With a further increase in the reaction rate or the sulfite concentration, respectively, most of the oxygen is already consumed close to the gas-liquid interface. Note that this observation is very similar to the results of Kück et al. [14], where for a constant sulfite concentration the reaction rate was modified via the water temperature.
To determine the mass transfer coefficient and also to have a better visual comparison to the numerical results, the gray-level images were evaluated based on the Stern-Volmer relationship. Figure 10 shows the resulting quantitative oxygen distributions. Note that concentrations below 0.1 mg/L are depicted white. This threshold also removes the areas with strong reflections close to the gas-liquid interface, and therefore, parts of the boundary layer where the oxygen concentration is presumably higher than the lower threshold. Additionally, the images are cropped to a size of 80 × 783 pixels.
The pseudo-color images in figure 10 reveal more details than their gray counterparts. Taking a close look at the sub-figures c) and d) leads to the following observation: there are two unexpected regions of higher concentration. The first one is close to the interface in proximity of the bubble south-pole. The second region is farther downstream before the recirculation zone ends and the wake becomes closed again. Based on the simulations, one would expect the highest concentration to be measured close to the ring where the flow detaches from the interface since the oxygen transferred into the boundary layer accumulates there. With increasing distance from the interface, the measured concentration of the transfer species (oxygen) should be always lower than close to the   interface since there is no mechanism to build up new concentration maxima 3 . Molecular diffusion or chemical reactions 4 can only lead to a depletion of oxygen. A possible explanation for the observed experimental concentration distribution will be given hereafter. As mentioned in section 3.3.2, we assume that the finite width of the laser sheet influences the measured concentration field. Figure 11 shows slices through the numerically obtained concentration fields extracted at different distances from the bubble center. Looking at the first column, the oxygen distribution is very similar to the experimental result, reflecting the ring-vortex. Because of the immobility of the interface, the contact time between fluid elements and the interface is significantly increased. Even for low Damköhler numbers like Da = 0.1, the oxygen in the wake is quickly consumed (the characteristic reaction timescale is 10 times larger than the hydrodynamic timescale). For the last column in figure 11, it is much harder to comprehend the concentration field. What we see is more like a slice through the boundary layer than through the wake.  Figure 11. Slices extracted at different distances z (in mm) from the bubble center (numerical data). Figure 12 shows the oxygen distribution in the wake obtained if the post-processing procedure described in section 3.3.2 is applied. The images in each row from left to right present the averaged concentration field for an increasing width of the laser sheet (averaging volume). For example, images in the first column "two layers" correspond to an average over the volume between the first and second slice in figure 11, or to a laser sheet thickness of 2 · 0.05 mm = 0.1 mm (the factor 2 results because we assume symmetry with respect to the bubble center). The image in the last column corresponds to the laser sheet thickness in the experiment. Interestingly, the averaging leads to a region of higher concentration at the end of the recirculation zone as observed in the experiments. Another effect is that the oxygen field close to the interface is blurred, but this area is usually not visible in LIF experiments because of reflections. Possibly one could find even better agreement with the experimental results by varying the averaging volume. For example, averaging only over the volume between the last three or four columns in figure 11 would presumably lead to a region of high concentration close to the bubble south pole. Further aspects which we have not considered here are that (i) in an experiment it is highly unlikely that the laser sheet insects the bubble perfectly centered, and that (ii) the bubble path always undergoes small oscillations. These effects could further influence the observed results.

Local and global mass transfer.
To determine the mass transfer coefficient k L , we follow the procedure of Kück et al. [13,29,14], and Jimenez et al. [11], but we would like to provide a more detailed derivation which clarifies assumptions and limitations. The global liquid-side mass transfer coefficient k L is defined as the molar fluxṄ Σ over the gas-liquid interface Σ per transfer area A Σ and macroscopic concentration difference ∆c O 2 , i.e., The concentration difference ∆c O 2 = c * O 2 −c O 2 ,∞ is defined as the difference between the saturation concentration of oxygen at the liquid side of the interface and the bulk concentration far away from the bubble (here we assume c O 2 ,∞ = 0 mol/m 3 ). Assuming a steady, spherical bubble shape, the effective mass transfer area is A Σ = πd b 2 . The only missing quantity to compute k L is the molar oxygen fluxṄ Σ . If the concentration gradient at Σ was known, the oxygen flux could be computed directly via Fick's law. Unfortunately, it is not possible to evaluate local concentration gradients  The temporal evolution of c O 2 (t, x) within Ω l (t) may be computed using the Reynolds transport theorem: where ∂Ω l (t) is the boundary of Ω l (t), u b (t, x) is the velocity of ∂Ω l (t), and n b is an outwardpointing unit vector normal to ∂Ω l (t). The concentration within Ω l (t) is time-dependent and may change due to convection, diffusion, and chemical reactions. Therefore, the partial time derivative is to be replaced by where ∇ denotes the Nabla operator, and u(t, x) is the velocity vector in the liquid phase. Substituting equation (11) in (10) and applying the divergence theorem leads to If we consider a co-moving domain of fixed shape Ω l and assume all processes to be steady, the temporal change of c O 2 within Ω becomes zero. Furthermore, the integral over the boundary ∂Ω l may be decomposed into two parts: the gas-liquid interface ∂Ω Σ and the outer boundary ∂Ω o . Simplifying and reordering equation (12) leads to If Ω l is co-moving with the bubble's center of mass, the velocity normal to ∂Ω l is equal to the bubble rise velocity u b · n b = U b on ∂Ω Σ . On the outer domain boundary ∂Ω o , the sign changes because normal and velocity vector point in the opposite direction: u b · n b = −U b . Because the bubble shape is steady, the interface velocity must be equal to the bubble rise velocity such that (u b − u) · n b = 0 on ∂Ω Σ . If the outer boundary if far from the bubble, the liquid velocity is at rest again such that u = 0 on ∂Ω o . The integral over the diffusive flux normal to ∂Ω Σ is the molar fluẋ N Σ required to compute the mass transfer coefficient. The diffusive flux over the outer boundary may be assumed to be negligible since the wake is closed and the oxygen concentration changes only mildly in streamwise direction; see figure 10. After applying all aforementioned simplifications, equation (13) becomesṄ Note that the negative sign of the oxygen flux means that it enters the control volume Ω l because of the outward-pointing normal vector. Expression (14) states that molar flux transferred from the gas to the liquid phase is equal to the sum of oxygen consumption within Ω l due to chemical reactions and the oxygen leaving Ω l over the outer boundary ∂Ω o due to the motion of the boundary. The planar LIF measurement provides concentration values within a plane through Ω l , which can be used to compute the integrals in equation (14) under the assumption of axis-symmetry around the rise velocity vector of the bubble. However, the computation of the volume integral over the reactive term in (14) poses new challenges compared to previous works: a significant amount of oxygen will be consumed within the boundary layer, but the boundary layer cannot be captured accurately. Thus, it is unlikely to obtain meaningful reactive mass transfer coefficients based on the present experimental approach. For experiments without chemical reaction, the surface integral can be computed as demonstrated in [14]; see figure 5 in the reference for a detailed sketch.
Since the experimentally observable area is limited, it is of interest to know at what distance behind the bubble the perturbation of the velocity field has vanished. Therefore, we extracted the velocity component in rise-direction on the centerline in the bubble wake (numerical data). Figure  14 shows the normalized velocity component plotted against the normalized distance from the bubble south-pole. At distances of y Σ = 2.45d b for the clean and y Σ = 3.01d b for the contaminated cases, the velocity field is perturbed by less than 1 % relative to the rise velocity. The Marangoni forces introduce additional shear forces into the liquid bulk, which perturb the liquid bulk more than in the clean case. The length of the wake perturbation scales approximately linear with the bubble diameter for clean and contaminated cases.  Figure 14. Ratio of the absolute velocity component in rise direction |u y | normalized with the bubble rise velocity U b plotted against the dimensionless distance to the gas-liquid interface y Σ /d b on a center line (in rise direction, y) in the bubble wake. The distances y Σ (u b = 0.99U b,clean ) = 2.45d b and y Σ (u b = 0.99U b,cont ) = 3.01d b depict at which distance the liquid phase is close to the steady rise velocity in the moving reference frame or close to be at rest in a fixed frame of reference, respectively. Figure 15 depicts, similar to figure 14, the dependency of the measured mass transfer coefficient on the distance from the bubble south-pole. After a distance of approximately 1.5d b , the k L signal becomes stable. Note that the mass transfer coefficient in the considered diameter range should be almost constant according to the literature references and the simulations. The dashed line in figure 15 depicts an exemplary reference. However, the computed k L values are not monotonously changing with the diameter, and the deviation from the reference is much larger than the confidence interval. The standard deviation in the figures 15 and 16 only reflects the changes in the measured concentration profiles. It is likely that the offset stems from errors in the measured diameter and rise velocity. The standard deviation of the rise velocity is larger than that of the diameter (see table 3), but the error in d b influences the computed k L quadratically; see equation (9). Furthermore, there is probably a systematic error in the measured concentration profile due to the finite thickness of the laser sheet, as discussed in section 4.2.  Figure 15. Computed mass transfer coefficient k L plotted against the normalized distance from the interface. The moving average and standard deviation (colored lines) were computed on a window size of 20 samples. The shaded area around the moving average expresses the 68% confidence interval (k L ± σ k L ).
For the k L values depicted in figure 16, the last 250 lateral planes were averaged over the measured concentration profiles in the bubble wake. Because the mass transfer coefficient in the considered diameter range is almost constant, we also computed an average value with respect to the diameter (a linear trend line would yield a similar result). Despite the different error sources, the diameter-averaged mass transfer coefficient agrees very well with Frössling's correlation and also with the numerical predictions. The relative standard deviation of approximately 20 % seems reasonable given the complexity of the overall measurement; see table 4 for the exact values. It is hard to reliably compare this value to other works because the overall confidence interval is usually not reported. Frössling mean k L exp.
68% CI simulations exp. Alves experiment Figure 16. Mass transfer coefficient from experiments, simulations, and literature plotted against the bubble diameter. The reference line is computed based on the drag coefficient by Schiller and Naumann (table 7, solid particle) and the Sherwood number correlation by Frössiling (equation (2)). Alvis experimental results are extracted from figure 10 in [1]. Table 4. Minimum, maximum and mean relative standard deviation for the computed k L .
In table 5 we also reported the numerical results for mass transfer at a clean bubble. A comparison with experimental data is not possible due to the inevitable presence of the fluorescence tracer. However, the deviation from different literature correlations is marginal. A modification of the coefficient in equation (18) fits the numerical results even better: Note that this modification is not meant as a "correction" but rather a step to avoid a propagation of the deviation into the reference curves for the reactive Sherwood numbers reported later on in figures 17 and 18. The numerically computed mass transfer coefficients for the contaminated cases agree very well with Frössling's correlation. It is remarkable that the sheer influence of Marangoni forces on the boundary layer structure leads to a more than 6 times lower mass transfer coefficients/Sherwood numbers. A minor modification of equation (2) that better fits the numerical prediction reads Sh = 2 + 0.5847Sc As we reported in section 4.2, already a slow chemical reaction leads to a substantial depletion of the oxygen concentration in the bubble wake. Because of the reason mentioned above, it is Table 5. Predicted k L /Sh for a mobile interface: present simulations, correlation by Takemura and Yabe [23] and correlation by Lochiel and Calderbank [16]. Both correlations can be found in A.2. The Reynolds number to evaluate Sh(Re, Sc) was computed based on the numerically predicted rise velocity.  (table 7) and the modified version of Lochiel's correlation for mass transfer (15). not possible to measure meaningful reactive mass transfer coefficients. However, we can provide insightful numerical results which we compare to the enhancement predicted by film-theory (3).
For the fully mobile interface the increase of the Sherwood number is small, E = Sh r /Sh ≤ 1.1; see figure 17 for the investigated reactive time scales. In this range, the relative difference between film theory and simulation is less than 1 %.
For the partially immobilized interface the reactive enhancement is much stronger E ≤ 2.6; see figure 18. We investigated the same ratio of convective to reactive timescale Da, meaning that we decreased the reaction rate corresponding to the decrease in the rise velocity. This adjustment is necessary to isolate the effect of surface contamination on mass transfer. In an experiment comparing clean and contaminated bubbles one would probably keep the reaction rate constant, such that Da increases significantly for the slower rising, contaminated bubbles. The observed reactive enhancement would be therefore even stronger than reported here.
The film theory predictions deviate up to 6 % from the numerical results. An empirical modification of equation (3) that fits the numerical results within a 1 % range reads  (table 7) and a modified version of Frössling's correlation for mass transfer (16). The fitting corresponds to equation (17).
Finally, we would like to explain why the film theory prediction of the enhancement factor for clean bubbles is more accurate than for contaminated ones. Figure 19 shows the local Sherwood number plotted over the polar angle. As already observed for the integral Sh values, the enhancement at the immobilized interface is much stronger. The main difference, however, results from the non-uniform local mass transfer enhancement in the presence of surfactant. Especially in the rear part of the bubble, where the contact time between fluid-elements and interface is the largest, the enhancement is stronger than at the upper bubble half. Agreement with the film theory can only be achieved when the enhancement is uniform and small, like for the fully mobile interface.  [13,29,14], but the source could not be isolated. Jimenz et al. [11] investigated the impact of (intentionally added) surfactants on mass transfer at rising bubbles. The authors of [11] also measured the influence of the same fluorescent dye (by the same provider) as used in this work on several system properties; see table 1 on page 113 in the reference. While there is no measurable influence on the water's density or viscosity at room temperature, the surface tension was reduced by approximately 0.5 %. This marginal change seems to be negligible at first sight. However, the difference of the flow around a rising bubble and the reduction of the rise velocity in a contaminated system mostly stem from Marangoni-forces, which are equal to the surface tension gradient and not to the surface tension itself. Simply said, a small change of surface tension that occurs over a small surface region will have a similar effect as a larger reduction happening over a larger region. Thus, in the surface tension dominated regime, a marginal change of the surface tension may result in a substantial difference of the macroscopic behavior. In fact, the difference between the maximum and minimum surface tension in the numerical simulation is only about 0.9 %. This change of the jump-conditions at the interface leads to a rise velocity reduction of about 50 %. Despite the evidence presented before, we cannot assure that the fluorescent dye is the surface active agent. There are several issues which prevent us from drawing clear conclusions: • Tap water was used in the experiment. Because tap water naturally contains contaminants, there could be the effect of several surfactants acting at the same time. The difference in the recorded rise velocities with and without fluorescent tracer could also stem from different degrees of contamination in the tap water since the experiments were not conducted on the same day. • The surface tension reduction, measured by Jimenz et al. [11], is very small and could have been below the accuracy range of the measurement device. • The purity of the fluorescent dye is only 98 %. The additional impurities could alter the surface properties, too 5 . Possible impurities are ruthenium chloride (RuCl 3 ) and biby (2,2-Bipyridine 6 ) which are left over from the synthesis of the fluorescent tracer. In water, the ruthenium chloride will dissociate into ruthenium and chloride ions (Ru 2+ , Cl − ), which are not surface active. The bipy, however, is hydrophobic and could adsorb onto the interface. Despite the ambiguous arguments presented before, it is not unlikely that the fluorescent dye caused the drastic drop of the rise velocity in the experiments. Therefore we discuss its molecular structure briefly hereafter.  [28]. Figure 20 depicts the spacial structure of the fluorescent dye (crystalline structure). The central ruthenium atom (turquoise) is twice positively charged. This charge is distributed by the organic phenanthroline ligands surrounding the central atom (blue -nitrogen, dark gray -carbon, light gray -hydrogen). Somewhere in the space around this molecule, there are two negatively charged chloride ions (not depicted) to complete the charge balance. Both the ruthenium complex cation as well the chloride ions are charged, but the charge density of the chloride is much larger since the atoms are much smaller. Moreover, the organic ligand has hydrophobic properties. Compared to the "classic" structure of surfactants consisting of a "polar" head (hydrophilic) and a non-polar "tail" (hydrophobic), the fluorescence tracer could behave like a hydrophobic sphere (the ruthenium complex cation) with one or two polar heads (the chloride ions). Thus, it could potentially be surface active. To draw definite conclusions, further experiments under well-defined conditions are necessary. A very sensitive way to investigate the adsorption of a chemical species onto the interface is the measurement of the bubble rise velocity over time (so-called local velocity profiles [19]). Under the influence of soluble surfactant, the rise velocity will show three very distinct stages of acceleration, deceleration, and steady state [20]. Of course, such experiments would have to be carried out using ultra-purified water and a fluorescent tracer with higher purity.

Conclusion
In this work, we investigated the reactive mass transfer from small rising oxygen bubbles by means of experiments and numerical simulations. The high-speed planar-LIF technique is an excellent tool to determine the mass transfer coefficient of such small bubbles. Based on the averaging procedure of several ROI-images, the signal to noise ratio could be improved significantly, such that a more accurate determination in comparison with [13,29,14] was possible. The fluorescent dye used in the experiments led to a strong reduction of the observed rise-velocity as it is typical for surfactants. We therefore discussed the potential polarity of the ruthenium complex and included the transport and surface tension reduction due to a surfactant in the numerical method. The simulations proved to have reached a high level of accuracy which allows to determine the bubble hydrodynamics and mass transfer for realistic time-and length-scales. Thus, simulations can validate experimental results and vice-versa. They also provide valuable insights into local quantities close to the interface.
There are several possible improvements to increase the quality of the experimental and numerical techniques further. In the experiments, a reduction of the uncertainty in each quantity involved in the calculation (diameter, rise velocity, wake concentration) will ultimately lead to a narrower confidence interval for the mass transfer coefficient. Future work will also be focused on employing novel reaction systems, such that educts and products can be detected simultaneously. On the simulation side, the subgrid-scale modeling has to undergo further developments to handle coupled higher-order reactions. Also, a mesh-topology with more control-volumes concentrated in the bubble wake would yield sharper concentration fields at larger distances from the interface. compared to bubbles rising in purified (deionized) water (see for example figure 7.3 on page 172 in [5]).
Takemura and Yabe [23] modified relation (18) based on their experimental data in order to cover a broader range of Reynolds numbers. Note that relation (22) in the reference is miss-printed. The correct correlation reads A. 3. Reactive subgrid-scale model. For this work, we extended the non-reactive subgrid-scale model [20] to handle convection-dominated species transfer accompanied by a first-order decay reaction. The non-reactive model was first introduced in [31] and then included in the interfacetracking framework [20]. A first step towards a reactive model was performed by Gründing et al. [8]. Here we only outline the main steps of the extension. A detailed description of the basic ideas and algorithms can be found in the cited literature. The substitute setting sketched in figure 21 for the decay reaction is the same as for the non-reactive model. The initial boundary value problem describing the reactive boundary layer consists of the partial differential equation accompanied by the initial and boundary conditions c(x > 0, y = 0) = 0; c(x = 0, y) = c Σ ; c(x → ∞, y) = 0 .
with the abbreviations In equation (22), erfc denotes the complementary error function. The next step in the modeling process is to adjust the analytical solution to the cell-centered concentrationc num = c num /c Σ of each interface cell in the numerical simulation. More specifically, the condition to compute the model parameter y reads where L is the thickness of the interface cell. An efficient iterative algorithm to find y is a Newton method based on equation (24). The following two expressions are necessary for the implementation: After y has been computed, convective and diffusive fluxes in the interface cells can be corrected. Equation (22) can be used directly to get a more accurate approximation of the concentration value at x = L. For the diffusive fluxes at x = 0 and x = L two more relations have to be used: 23 The chemical reaction introduces another difficulty compared to the non-reactive model. As y becomes larger, the boundary layer and its thickness approach a steady state. In the limit y → ∞ the solution of the simpler film-theory is obtained: Consequently, also the average model concentration value in relation (24) is limited bȳ If the cell-centered concentration valuec num is larger than expression (30), it becomes impossible to find a suitable y. Such values typically appear in the rear-part of a rising bubble. The normal velocity component close to the interface causes an additional convective transport of the chemical species, and hence the boundary layer becomes thicker. Because normal velocity contributions are not included in the substitute problem, the model function cannot be used without modification. A heuristic solution which works excellently consists of the following two steps: (1) When the maximum average concentrationc max is smaller thanc num we fit the film theory solution by adjusting the parameter A. Decreasing A corresponds to increasing the ratio of diffusive flux into the boundary layer to reactive decay. This adjustment emulates the additional convective flux normal to the interface in the numerical simulation. (2) We do not apply the full SGS-model fluxes but rather a combination between standard (linearly interpolated) and SGS-model fluxes. A simple blending of the form φ = (1 − w)φ SGS + wφ num works well with φ representing either the concentration value or its gradient. A sensible weight is w =c num . As c num approaches the interface value c| Σ , the boundary layer is typically well resolved and no model correction is required. On the other hand, when c num is very small, a simple linear approximation will overestimate fluxes leaving the interface cell, while the SGS-model fluxes yield much better results. A standard test case for SGS-model validation is the semi-analytical solution of species transfer from a spherical bubble rising at a very low Reynolds number. The solution is termed semianalytical because the species transport equation is solved numerically using an analytical expression for the velocity field to compute convective fluxes. Our main quality indicator is the local Sherwood number depicted over the polar angle. Figure 22 highlights two points: (i) The species transfer enhancement due to the chemical reaction in the boundary layer is captured well, and (ii) the solution shows very little dependence on the computational mesh. The reference solution was computed on a mesh where the first cell layer thickness was approximately 0.1 µm. Even the global Sherwood numbers computed on the coarsest mesh, employing the subgrid-scale model, deviate by far less than 1 % from the reference. We also repeated the validation procedure for P e numbers up to P e = 10 6 , but the results were similar and are therefore not included here.