Numerical and experimental analysis of oxygen transfer in bubble columns: Assessment of predicting the oxygen-transfer rate in clean water and with surfactant solutions

The purpose of this study was to develop a numerical model to estimate the oxygen-transfer rate for a laboratory-scale bottom aeration system at a 1.28 L reactor volume and to contribute to fundamental knowledge regarding the oxygenation of surfactant solutions. The primary goal of the study has been to develop a computational fluid dynamics (CFD) model using Euler – Euler (EE) mixture model coupled with the advection-diffusion equation to predict the oxygen-transfer rate in bubble columns containing clean water. The secondary goal has been to apply the model to water-based solutions containing the surfactant lauric acid (DDA) to identify options for further development of the model to make it applicable to surfactant solution systems. The Sauter mean diameter (SMD) was calculated to represent the average bubble diameter, based on available experimental data for different combinations of superficial velocities rate and DDA concentration. The oxygen-transfer rate in clean water fit well with experimental data at lower superficial velocities, and the differences in volumetric mass-transfer coefficients were 0.7% and 3.3% for superficial velocities of 0.24 cm/s and 0.48 cm/s, respectively. For surfactant solutions, the model overestimates the oxygen-transfer rate due to surfactant adsorption at the bubble/water interface and the consequent decrease in the mass-transfer coefficient not being modeled sufficiently. A correction factor for the mass-transfer coefficient based on a larger sample size of experimental data may need to be calculated and applied to improve model predictability.


Introduction
Contaminants vary between municipal wastewater and industrial wastewater, as well as between different industrial effluents, and so the treatment techniques required vary accordingly.The pollutants in municipal wastewater are mainly relatively easily degradable organic compounds, which are typically removed by aerobic biological treatment techniques, mainly the activated sludge process (Kehrein et al., 2020).The pulp and paper industry produces wastewater with similar biodegradable organic matter, and aerobic biological treatment techniques are most often used for its treatment.Typically, variants of the activated sludge process or aerated ponds are used (Ashrafi et al., 2015).However, biological purification methods are not suitable for all types of pollution; wastewater treatment plants, both municipal and industrial, usually include several different treatment steps to deal with all the different substances that may be encountered.Aerobic biological purification techniques are based on aerobic bacteria, which rely on oxygen as an electron acceptor when breaking down organic pollutants in water.The rate of degradation and efficiency of the process are closely linked to the concentration of dissolved oxygen (DO), which is one of the central operating parameters analyzed in controlling the purification efficiency of the outflow from aerobic biological purification stages (Garcia-Ochoa and Gomez, 2009;Pittoors et al., 2014aPittoors et al., , 2014b)).In order for oxygen to dissolve in wastewater at a sufficiently high rate to meet the oxygen demand of bacteria and thereby maintain a high rate of decomposition of pollutants, aeration techniques are required.Aeration further contributes to the mixing of wastewater and sludge (the biomass of bacteria), which in turn increases the purification efficiency (Baquero-Rodriguez et al., 2018).The most common aeration technique is bottom aeration, in which air is pumped through pipes to diffusers placed near the bottom of the purification reactor, from which it emerges through small openings to produce bubbles with high specific boundary area, by which oxygen can be transported from the air to the wastewater (Pittoors et al., 2014a(Pittoors et al., , 2014b;;Amaral et al., 2017;Drewnowski et al., 2019).

Nomenclature
Active sludge technologies have been developed and optimized over a long time, and have consistently demonstrated high purification efficiency for organic matter and nutrients, but an associated disadvantage is the high energy demand that aeration entails.Activated sludge processes have been reported to account for up to 1% of the total electricity consumption in developed countries, and bottom aeration is by far the most energy-intensive aspect (Meerburg et al., 2016).Indeed, bottom aeration can account for 45-90% of the total electricity consumption of a treatment plant, depending on the size of the operation and the technique used (Baquero-Rodriguez et al., 2018;Amaral et al., 2017;Drewnowski et al., 2019;Karpinska and Bridgeman, 2016).Accordingly, bottom aeration can account for up to 49% of the total operational cost (Drewnowski et al., 2019).
Extensive knowledge of the factors that influence parameters such as the oxygen-transfer rate (OTR) from air bubbles to wastewater; the oxygen uptake rate (OUR); oxygen-transfer efficiency (OTE), that is, the proportion of the oxygen that leaves the diffusers that has time to dissolve in the water, and the DO content, is important in order to optimize the aeration process and minimize energy consumption (Pittoors et al., 2014a).For example, continuous monitoring of the DO and appropriate adaptation of the airflow has been shown to reduce electricity consumption for aeration by 26% (Guven et al., 2019).
A common class of pollutants in both municipal wastewater and wastewater from the pulp and paper industry comprises surfactants.Surfactants assemble to form interfaces between substances, as well as interfaces between thermodynamic phases of the same substance (Palmer and Hatley, 2018).Surfactants are present in many everyday products, such as soap and detergents, and after household use eventually end up in municipal treatment plants, and also take the form of fatty acids in wastewater from the paper and pulp industry (Ashrafi et al., 2015;Baquero-Rodriguez et al., 2018;Jimenez et al., 2014).The surfactant concentration in municipal wastewater is typically in the range 10-20 mg/L, but can reach 300 mg/L in some industrial wastewaters (Dereszewska et al., 2015).The presence of surfactants has been shown to negatively affect oxygen transport when aerating wastewater, as a result of their assembly at the interface between air bubbles and water restricting oxygen transport (Chen et al., 2013;Rosso and M, 2006).Contrasting results have been obtained on how the volume fraction of air in bubble columns is affected by surfactants, with some studies indicating that the volume fraction increases with the surfactant concentration and others indicating no such relationship.Prakash et al. (2021) observed an increasing volume fraction of air in a reactor at all tested airflows, which they attributed to reduced surface tension with increased surfactant concentration leading to the production of smaller bubbles at a lower rate.However, Belo et al. (2011) detected no difference in the volume fraction of air with increasing surfactant concentration, although a decrease in bubble diameter was identified.
In all relevant studies, the volumetric mass-transport coefficient (k L a) has been shown to decrease with increasing surfactant concentration (Chen et al., 2013;Painmanakul et al., 2005;Belo et al., 2011;Painmanakul et al., 2005;Vasconcelos et al., 2003).Several studies have shown that the decrease in k L a is not further affected by increasing the surfactant concentration above a certain threshold (Belo et al., 2011;Vasconcelos et al., 2003).Chen et al (Chen et al., 2013).observed that k L a recovered slightly at higher surfactant concentrations as the specific boundary surface area (a) continued to increase slightly while the mass-transport coefficient (k L ) stabilized.Vasconcelos et al (Vasconcelos et al., 2003).noted a flattening off of the k L a evolution plot with increasing surfactant concentration beyond the critical micelle concentration (CMC).In several studies, it has been shown that a increases with increasing surfactant concentration (Chen et al., 2013;Painmanakul et al., 2005;Belo et al., 2011;Sardeing et al., 2006).This follows from the fact that both the bubble size and the bubble rise velocity decrease with increasing surfactant.Belo et al. (2011) pointed out that k L must diminish at a faster rate than the increase in a in order for k L a to decrease with increasing surfactant concentration.Both Chen et al. (2013) and Belo et al (Painmanakul et al., 2005).also observed a decrease and levelling off the rate at which a increases with increasing surfactant concentration.Studies have consistently shown that k L decreases with increasing surfactant concentration (Jimenez et al., 2014;Chen et al., 2013;Belo et al., 2011;Painmanakul et al., 2005;Vasconcelos et al., 2003;Fletcher et al., 2017), eventually levelling off to a constant value (Belo et al., 2011;Vasconcelos et al., 2003).Vasconcelos et al. (2003) showed that the value of k L consistently falls between the theoretical value according to Higbie's penetration theory for mobile interfaces and Frossling's theory for rigid interfaces.In contrast to the aforementioned studies, the system analyzed by (Gomez-Diaz et al., 2009) showed an initial decrease in k L with increasing surfactant concentration, followed by an increase.
The oxygen-transport rate during bottom aeration is influenced by a large number of factors besides the presence of surfactants, including geometric parameters (e.g., reactor dimensions, diffuser depth, and bubble size), diffuser type, physicochemical parameters (wastewater chemical composition, temperature, viscosity, etc.), and flow-related parameters (Pittoors et al., 2014a(Pittoors et al., , 2014b;;Baquero-Rodriguez et al., 2018;Al-Ahmady, 2011).The relationships between these influencing parameters and the oxygen-transport rate are complex, making it difficult to develop universal empirical correlation models to predict oxygen-transport rates in all systems.The use of numerical methods to deal with this complexity in systems ranging from simple laboratory-scale reactors to full-scale active-sludge pools has increased rapidly in recent decades.Such systems are characterized by a variety of physical, chemical, and biological phenomena occurring simultaneously.As a result, several numerical models are often linked together, such as CFD for flow, advection-diffusion models for oxygen transport in water, and the activated sludge model (ASM) for biological activity (Karpinska and Bridgeman, 2016).
Due to wide-ranging flow behaviour in bubble columns, the structure of hydrodynamic flow regimes is an important to understand in the design and scale-up of bubble column reactors (Clift et al., 1978).A successful way to do this is to combine experimental work with numerical simulations (Shaikh and Al-Dahhan, 2007).A multitude of numerical simulation studies have been performed on bubble columns, encompassing fluid flow and mass transport, with great variation in the models for turbulence, the two-phase system, and forces acting on the interfaces between the phases.However, very few studies have been carried out with extensive assessment of the effects of surfactants on hydrodynamics or mass transport.Some CFD studies have been conducted on bubble columns with water/surfactant solutions, which differed from CFD studies on bubble columns with pure water in that an empirically derived constant was invoked in the equation to correct for resistance, thereby giving better consistency with experimental results for the air volume fraction in the reactor (Ertekin et al., 2021;McClure et al., 2015).Ertekin et al (Al-Ahmady, 2011).and Fletcher et al (Fletcher et al., 2017;Ertekin et al., 2021) used empirically derived values for k L in studying the effects of surfactants on oxygen transport.However, the authors of all three aforementioned studies agree that more research is needed to develop customized equations and correlations to address the effect of surfactants on resilience and k L that are applicable to a wider range of systems.Bastani et al. (2018) studied the dynamic surface tension (and surface tension gradient) of a single bubble by CFD, invoking relationships for the transport of surfactants from water to the interface with a single bubble, and along the interface, as well as relationships for surfactant adsorption and desorption at the surface.The authors highlighted the value of the CFD study in providing insight into the local concentration profile along the surface, which is not possible from an experimental standpoint.
To guarantee sufficient purification efficiency in an aerobic biological purification system as energy efficiently as possible, knowledge of the oxygen transport in the system is required.When defining new systems and calibrating existing systems, tools are needed to predict oxygen transport so that a working solution can be identified and applied.The purpose of this study is to facilitate the prediction of oxygen-transport rate during bottom aeration in laboratory-scale systems, and to contribute to the knowledge base regarding oxygenation of aqueous surfactant solutions, more specifically those containing lauric acid (DDA).Thus, a goal of this study has been to develop a numerical model with short simulation times to recreate the oxygenation process in a bubble column reactor on a laboratory scale.The primary sub-goal has been to simulate the oxygenation process and compare it with experimental data to obtain a model that can predict the oxygen-transfer rate in the case of pure water.The secondary sub-goal has been to use the model to simulate the oxygenation process for water containing different concentrations of the surfactant DDA, to analyze the differences between simulated results and experimental data, and to identify how the model needs to be further developed in order to be able to characterize the oxygenation process in surfactant solutions.

Materials and methods
The method section is divided into a first part concerning the numerical modeling, including the geometry and mesh of the model, as well as the use of mathematical relationships, assumptions, boundary conditions, and initial states.In the second part, the procedure for bubble measurements and the ensuing calculations is presented.
Data were obtained from experiments performed in the Department of Chemical and Engineering Sciences, Karlstad University.These data were used for validation and in some cases for input to the model.Data were obtained in the form of the oxygen concentration versus time during the aeration experiments and volume fraction of air in the reactor.Moreover, photographs were taken of the reactor during operation.Experimental oxygen concentration data were used solely for comparison with the results of the simulation: a direct comparison and to calculate and compare k L a. Data on the volume fraction of air were used for comparison and to adjust the reactor height in the model.The volume of air was measured by considering the additional reactor volume when the aeration approached steady-state conditions (Fig. 1).The photographs of the reactor were used to estimate the bubble sizes, which underpinned the bubble diameter equations used in the model.Data were obtained from experiments with airflow rates of 0.1, 0.2, and 0.3 L/min (equal to superficial velocities at 0.24, 0.47 and 0.71 cm/s) in pure water and in DDA solutions in the concentration range 0.3125-15.625mg/L.Lauric acid, with the IUPAC name dodecanoic acid (DDA), is found in forest industrial effluents, and was therefore selected as a model substance.The relevant experiments were performed in a reactor of diameter 30 mm and containing 1.28 L of water, hence resulting in a reactor height at approximately 1.846 m.

Modeling and simulation process
The bubble column model was built in COMSOL Multiphysics 5.5, specifically the modules "Bubbly Flow", which is an Eulerian-Eulerian (EE) mixture model, and "Transport of Diluted Species", which is based on an advection-diffusion equation.The simulations were time-dependent and a duration of 500 s was simulated for all combinations of superficial velocities and DDA concentration.The flow was turbulent and URANS equations with standard k-ε were used as the turbulence closure model.The default values of the module were retained as the coefficients in the turbulence equations.Simulations were performed for three superficial velocities in pure water and for an additional 23 combinations of the same superficial velocities with varying DDA concentrations.Selected combinations were simulated to enable comparison with available experimental data.

Eulerian-Eulerian mixture model
The dynamics of the bubble and water flow was solved by the EE mixture model, which is a turbulent two-phase model.This physical approach solves the single momentum of the mixture, the continuity equation for the liquid and gas phases, and the gastransport equation [Eqs.( 1)-( 3)].
∂ ∂t The gas side velocity u g,i is determined by the slip velocity term, u slip,i which gives the relative velocity between the gas and liquid phases.The term u drift,i , determined by the turbulent viscosity, is added to the gas-transport equation.
The drift velocity adds a diffusion term to the gas-transport equation.The term is implemented to reduce artificial accumulation of the gas volume fraction, which can appear in areas where the volume fraction of bubbles becomes stagnant.The increased diffusion of volume fraction provides a more robust numerical solution.Inserting Eq. ( 4) and Eq. ( 5) into Eq.( 3), the gas transport can be expressed by Eq. ( 6): The liquid and gas velocities are related by the pressure drag slip model, which is defined by Eq. ( 7) and is typically used to estimate friction drag on rising air bubbles (Chuang and Hibiki, 2017;Khan et al., 2020;Pourtousi et al., 2014;Shu et al., 2019;Sokolichin et al., 2004).
Here, the definition of the drag coefficient is determined by the formulation developed by Johansen and Boysen (Johansen and Boysan, 1988), Eq. ( 8): The Eö number (or Bond number) is the ratio of the gravitational to capillary forces according to Eq. ( 9): where the characteristic length is the average bubble diameter.Based on experimental trials, rough calculations of the Eö number suggested that capillary forces cannot be neglected.
The interfacial area is one of the parameters that determines the overall liquid-side mass-transfer rate.The physical presence of individual bubbles is not determined in the EE mixture model.Therefore, conservation of bubble number is set up according to Eq. ( 10): The number density of bubbles is specified at the inlet boundary, and the spatial distribution of these bubbles based on the gasphase velocity is calculated according to Eq. ( 10).
Bubble coalescence and splitting leads to both larger and smaller bubbles being formed relative to the average size in the column, which means that the bubble size distribution (BSD) becomes wider.Since most mathematical models that include bubble size, regardless of the phenomenon being analyzed, use a singular value rather than a distribution expression, a representative average value for BSD is required.The SMD (d 32 ) is used as the characteristic length scale to estimate the Reynolds number and is defined as the ratio of the volume equivalent diameter to the surface equivalent diameter according to Eq. ( 11) (Agrawal, 2013;Kowalczuk and Drzymala, 2016;Nekoeian et al., 2021).
where d v is the volume equivalent diameter.Rough calculations on the average bubble size and gas velocity from experiments suggest that the flow can be considered as fully turbulent.The standard two-equation k − ε model was used with an added bubble-induced source term according to Eq. ( 12): In the transport equation for the turbulent dissipation rate, the source term is expressed by Eq. ( 13) as follows: where C k and C ε are empirical constants.

The advection-diffusion model
The advection-diffusion equation is coupled to the EE mixture model.Hence, the species of interest (c), which is the transportation of the DO from the air bubbles to the water, is coupled with the velocity and pressure field solved in the EE mixture model.The transport of DO in water is expressed by Eq. ( 14): Here, the left-hand side of Eq. ( 14) represents the transient, diffusive, and convective term of the DO.The right-hand side is either a source or sink term, which in this system represents the mass transfer of oxygen from the gas side.The approach for modeling the source term is presented in the next section.

Oxygen source term
Bubbles are not represented physically in the computational domain, and therefore oxygen production has to be specified as a source term.The two-film theory, which uses Henry's law to compute the required mass transfer from gas to liquid, is expressed by Eq. ( 15): Here, c * s is the equilibrium concentration given by Henry's law according to Eq. ( 16): where p is the pressure computed in the EE mixture model, p 0 is the reference pressure, and H is Henry's constant, which is calculated according to Eq. ( 17).Here, h r is the height of the water column, y O2 is the molar fraction of oxygen in air, and M is the molecular weight of oxygen.The calculated value of H for oxygen in water was confirmed to be within the range reported in the literature (Sander, 2015).
The mass-transport coefficient k L is based on Higbie's penetration theory and is calculated according to Eq. ( 18), where D O2 is the diffusivity of oxygen in water, u slip is the relative velocity between air and water, calculated in COMSOL from the balance of pressure and resistance forces, and h b is the bubble height Nekoeian, Aghajani (Nekoeian et al., 2021), Hongprasith, Dolkittikul (Hongprasith et al., 2016), Kulkarni (Kulkarni, 2007), Higbie (Higbie, 1935).No data were found for the diffusivity of oxygen in DDA solutions, but since the solutions in the study were very dilute in terms of DDA, the diffusivity was assumed to be the same as that for pure water (1.97 ×10 − 9 m 2 /s).
Since bubbles are not physically represented in the CFD model, oxygen transfer is represented as a source term coupled to the interfacial area a, which is calculated in the Bubbly Flow module according to Eq. ( 19): The mass-transfer source term in Eq. ( 15) is accounted for in the advection-diffusion equations through the Reactions term: The magnitude of the mass transfer is dependent on the liquid-side mass-transfer coefficient and the interfacial area, which are determined from experimental trials.

Surface tension
The main physical effect of the presence of surfactants at interfaces is that they lower surface tension [17, (Nekoeian et al., 2021).To include the effect of surface tension variation with the DDA concentration, a polynomial regression model, Eq. ( 21), was adapted using the standard tool in Excel (R 2 = 0.9999) in order to calculate the surface tension as a function of the DDA concentration based on data published by (Kimura et al., 1998).The regression model is not based on any established physics, but was only developed to easily include the surface tension variations in the model of the DDA concentration range studied here.
The slip-velocity that is used in Eq. ( 7) is dependant of the drag form coefficient, which is dependent on the Eö number Eq. ( 9).The Eö number is explicitly dependant on the surface tension.The regression line is used so that when different DDA-concentrations are analysed, the drag coefficient will change which in turn will affect the slip velocity.

Initial and boundary conditions
The bubble column reactor was modeled according to the dimensions of the reactor used in the experimental trials.The geometry was built as two-dimensional axial symmetric, with a radius of 15 mm (i.e., a cylinder diameter of 30 mm) and a height at 1.846 m.The model did not include the initial movement upper end of the water volume (in reality, the total volume increases when aeration starts), so the reactor volume at various combinations of DDA concentration and superficial velocities was based on the corresponding experimental volume fractions of air in the reactor at maximum flow.
Boundary conditions were specified for the liquid phase and gas phase according to Fig. 1.Standard wall conditions were set on the solid surfaces, except for the free surface, for which a slip condition was set to enable tangential motion.A gas flux condition was specified at the inlet, which was determined from experimental trials and was ramped from zero to the set value to obtain a consistent initiation of the simulation.The air mass flow and the number of bubbles entering were determined according to Eqs. ( 22) and ( 23): Initial values for the simulation were set to zero vector for the water velocity field, zero for the volume fraction of air, 0.3 mg/L for the oxygen concentration in the reactor, 30 1/m 3 for the bubble density (the number of bubbles per m 3 ) to facilitate convergence, and COMSOL's default values for turbulent kinetic energy and energy dissipation.The initial pressure field was defined according to Eq. ( 24), where P ref is the reference pressure, ρ l is the density of water, g is gravitational acceleration, and z is the height coordinate.
Furthermore, a step function was defined to obtain a smooth function for the increase in air flux through the air diffuser (and the associated flux in terms of the number of bubbles) from zero to the fully developed value.The function was defined such that it reached its fully developed value after 4 s.Escalation was applied to avoid convergence problems resulting from large gradients at singular points.
For the liquid and air phases of the model, the materials predefined by COMSOL were "water, liquid" and "air", respectively, with associated physicochemical properties.The global temperature and reference pressure parameters were selected as 20 • C and 1 atm, respectively.The system was allowed to evolve isothermally and heat development from viscous energy dissipation was neglected.The DDA concentration in the study was low (dilute solution), and any changes in the density and viscosity of the solution could be neglected; the predefined material "water, liquid" was also used for DDA solutions.Transport of DDA to the bubble/water interface and the effects of dynamic surface tension were not included in the model.A gravitational node was added to include gravity in the transport equations.Three-dimensional simulations were performed at the superficial velocity 0.24 cm/s for all DDA-concentrations to test the assumption of rotational symmetry flow pattern yielded in the 2-dimensional axial symmetric domain.

Mesh-grid
The geometry was divided into two domains: one at the bottom (0 ≤ z < 0.05 m), where the velocity gradients were large, and the rest of the geometry (0.05 ≤ z ≤ hr m), where the gradients were smaller.In the bottom domain, an unstructured, triangular mesh calibrated for fluid dynamics was auto-generated.In the upper domain, a structured rectangular mesh was built, 60 elements in width and 180 elements in height.The width of the elements was kept constant, while their length was defined such that it increased linearly between two extremes, with the criterion that the length of the first elements in the positive z-direction (at z = 0.05 m) was 4% of the length of the last elements in the positive x-direction (at the free surface at the top of the reactor).The increase in the length of the elements of the structured mesh was included because the gradients were generally higher close to the bottom of the reactor.The mesh was also made finer at the boundary layer between the reactor wall and the water.

Bubble size measurement
Bubble size is a central parameter in bubble flow and has a major impact on their pressure and resistance balance (which affects the relative rate of transfer between bubbles and water), the flux in terms of the number of bubbles, the boundary surface area between bubbles and water, and the mass-transport coefficient.The model used in the study does not take into account any bubble size distribution.The bubble size is assumed to be uniform, which means that a good representation of bubble size is important.The parameter used in COMSOL to represent the bubble sizes is the bubble diameter.Here, the representative diameter used for the model was SMD (d 32 ), which was calculated based on experimental data.
Data concerning the shape and size of bubbles were collected from photographs of the bubble column reactor beside a tape measure during operation; see Fig. 2. To measure the size of the bubbles, images were imported into AutoCAD software, and the size was adjusted so that the ruler scale matched the software's defined length scale.Lines were then drawn along the larger and smaller axes of the bubbles (see Fig. 2), and all data on the lengths of the axes of the bubbles were exported and processed in Excel.Only those bubbles K. Rezk et al. that could be clearly discerned were included; overlapping bubbles and bubble swarms were excluded unless the contours of individual bubbles were clear.In total, 1865 bubbles were measured, distributed over 22 combinations of airflows and DDA concentrations.The curved reactor face was not considered during measurements.However, the distance from the picture taken in relation the position of the ruler is roughly 10 cm, which is the reason why it was neglected.
From observation, the bubbles usually take the form of ellipsoids; hence, a volume-equivalent diameter (the diameter of a sphere with the same volume as the bubble) was calculated as the representative diameter of the bubbles.The volume-equivalent diameter was calculated according to Eq. ( 25), where l short is the length of the shorter axis (m) and l long is the length of the longer axis (m).The depth axis was assumed to be of the same length as the shorter of the visible axes.The average representative diameter d 32 was then calculated according to Eq. ( 11), where n is the number of bubbles.
In an attempt to predict d 32 based on airflow and DDA concentration, an equation was developed that could be implemented in COMSOL to calculate d 32 .The calculated values of d 32 according to Eq. ( 11) followed a reasonably similar pattern for superficial velocities of 0.24 and 0.47 cm/s, whereas the d 32 of 0.47 cm/s deviated from the pattern.As there are large potential margins of error in the calculation of d 32 (e.g., uncertainties in the length of the invisible axis, bubbles without a clear contour, bubble swarms, and in some cases poor image quality), the values of d 32 for the superficial velocity of 0.47 cm/s were excluded and the analysis was based solely on the values of 0.1 and 0.3 L/min.No clear trend could be discerned for the values of d 32 at higher DDA concentrations.Due to the relative paucity of the data set for d 32 , a cubic spline interpolation was used instead of higher order regression models in order to avoid Runge's phenomenon, which is basically over-and under-shot of the polynomial curve interval between data points.Each cubic spline function was inserted in COMSOL so that the appropriate d 32 was chosen depending on the DDA concentration.However, more data at higher DDA levels, more specifically at 4 mg/L and above, are required to better understand how d 32 behaves.
To calculate the bubble height used for Higbie's penetration theory, the "Sauter mean" of the height (length in the vertical direction) was also calculated according to Eq. ( 26).The shorter axis of the bubbles was often almost parallel to the vertical direction, so the bubble height was assumed to be equal to the length of the shorter axis.The ratio of h 32 to d 32 showed no clear correlation to either superficial velocities or DDA concentration, so an average h 32 fraction of d 32 for the 22 combinations of superficial velocities and DDA concentration for which bubble photographs were available was used to calculate the bubble height according to Eq. ( 27).

Results and discussion
The simulated results for k L a for pure water at all superficial velocities are compared with experimental data in Table 1.The percentage differences for the three superficial velocities are 0.69%, 3.30%, and 14.03%, respectively.The numerical model uses Higbie's penetration theory to estimate the k L coefficient, which shows good agreement with the numerically estimated k L a values for superficial velocities of 0.24 and 0.47 cm/s.As the number of bubbles is increased, there is greater deviation of the homogeneous bubble distribution due to boundary surface blockage.This leads to bubble coalescence, which, in turn, limits the diffusive boundary layer.This explains the overestimation of k L a from simulated results at higher superficial velocities.
Unlike the results for k L a, the results from simulation for the air volume fraction in the reactor were best matched with the experimental data at the superficial velocity rate of 0.71 cm/s, whereas a difference of over 30% was obtained at the superficial velocity rate of 0.24 cm/s; see Table 1.Because the difference in k L a between the simulated results and experimental data was minor, as indicated in Fig. 3, experimental data for the volume fraction of air in the reactor at the 0.24 cm/s superficial velocities are assumed to be unreliable.The air volume fraction in the reactor is clearly correlated with k L a as the available specific surface area (a) depends on how many bubbles are in the reactor, and with a constant bubble diameter in the model, a 33% lower volume fraction of air would correspond to a 33% lower value of k L a, which was not the case here; see Fig. 3. Experimental data on the air volume fraction were collected by measuring the height difference in the reactor with a tape measure during active aeration, and the experimentally developed volume fraction of air in the reactor at the superficial velocity of 0.24 cm/s corresponded to a height difference of approximately 1.3 cm.This suggests large potential margins of error, as any measurement errors of a few mm can have a major impact on the result.
A comparison of DO concentrations over time between simulated and experimental data is presented in Fig. 3.Even though the k L a value is somewhat overestimated for a superficial velocity of 0.71 cm/s, the deviation of the aeration rate does not exceed 5%.
The results of the simulations and the experimental data for the relationship of k L a and DDA concentrations, and the corresponding fitted curves, are compared in Fig. 4. At all three superficial velocities, k L a decreases with increasing DDA concentration, except at very low concentrations, and the decrease appears to subside at higher concentrations, consistent with previously reported results (Belo et al., 2011;Vasconcelos et al., 2003).At a superficial velocity of 0.71 cm/s, the curve appears to flatten out at around 10-12 mg/L DDA, and even recovers slightly, as has been reported previously by (Chen et al., 2013).Experimental data show that even at surfactant concentrations below 10 mg/L, k L a is negatively affected by increasing DDA concentration, highlighting the importance of studying these effects in dilute surfactant solutions to optimize the purification technique, as surfactant concentrations in municipal wastewater are often in the range 10-20 mg/L, and can be much higher in industrial wastewater (Dereszewska et al., 2015).Comparison of the fitted curves for simulated and experimental data in Fig. 4 illustrates that the discrepancies between the k L a estimates increase with increasing DDA concentration.This further reiterates that due to the formation of larger bubbles, increased surface blocking, and non-uniformity, the calculated slip velocity ( u slip ) from the numerical model is not well represented in Higbie's penetration model, see Eq. ( 18).
The effects of DDA concentration in the numerical model were included through the expressions for d 32 [Eq.( 26)-( 28)] and the surface tension [Eq.( 21)].As shown in Fig. 2, d 32 clearly decreased at low DDA concentrations for superficial velocities of 0.24 and 0.71 cm/s and then started to increase, whereas the trend for d 32 at 0.47 cm/s showed minor fluctuations at low DDA levels.The decrease in d 32 and the following increase with increasing DDA concentration leads to the specific surface area (a) first increasing and then decreasing, which has a clear effect on k L a for the simulation.The variation of d 32 is clearly linked to the simulated results of k L a, see Fig. 4, in which the trend is the same but horizontally mirrored.
The volume fractions of air in the aqueous DDA solutions for the simulation and experimental runs are presented in Fig. 5. Good agreement is seen at low DDA concentrations, but the deviation increases with increasing DDA content and the simulated data overestimate the air volume fraction according to the trend lines drawn in Fig. 5.In the simulation, the volume fraction showed an increasing trend at higher DDA concentrations, whereas the experimental data showed a decreasing trend.Trends in air volume fraction at a superficial velocity of 0.24 cm/s were not analyzed further due to the inconsistent data representation in Table 1.The bubble velocity in the reactor determines the residence time of the bubbles, which directly determines the air volume fraction in a specific reactor with a defined superficial velocity.The variation in the volume fraction of air with the DDA concentration is hence directly affected by variations in bubble velocity.Two main factors influence the slip velocity between bubbles and water at varying   7)-( 9), the relative velocity decreases with decreasing surface tension for the entire DDA concentration range considered in this study, so the effect of surface tension is to increase the volume fraction of air in the reactor with increasing DDA concentration.
According to the theory established by Shu et al.( 2019), the flow resistance coefficient (C D ) of bubbles (leading to reduced relative velocity) increases with the addition of surfactants (Marangoni force), which is consistent with the effects of reduced surface tension on C D according to Eqs. ( 8) and ( 9).The results presented by Prakash et al. (2021) also showed an increasing air volume fraction in a bubble column when adding surfactants, which strengthens the simulated results, but contradicts the experimental data of this study.Since surfactants increase flow resistance and, in turn, decrease relative velocity, d 32 must be invoked to account for the increasing relative velocity with increased DDA concentration that results from the decrease in air volume fraction in the reactor.This needs to be further investigated experimentally by collecting more data for various reactor geometries to gain a better understanding of how the air volume fraction behaves at different DDA levels, given that it is critical to provide a better estimate of k L a.
The results on average u slip from the simulation are presented in Fig. 6 where an apparent reduction of the slip velocity on an average basis was observed for all superficial velocities and DDA concentrations.
The fully 3-dimensional computational domain was simulated in order to test the assumption of rotational symmetry flow pattern.As Fig. 7 is indicating, no major changes in are observed in the volumetric mass transfer coefficient, gas volume hold up and the average slip velocity.We believe that the reason is the diameter is small, hence, wall effect will keep the bubble distribution in a somewhat homogenous state.At larger diameters, the difference for these parameters should increase based on the axial-symmetric and fully 3-dimensional system.
Simulation results for DO, air volume fraction, slip velocity, and volumetric mass-transfer coefficient are presented in Fig. 8. Local heterogeneities of all parameters were observed at the bottom of the reactor.Moreover, a recirculation zone was observed at the lower end due to advective (higher water velocities) mass transfer followed by low radial dispersion of the air volume fraction, which led to reduced oxygen transport.
Since different surfactants work in different ways and with varying degrees of effectiveness, models need to start from more basic physicochemical properties of the surfactant solutions, including surface tension.Furthermore, the mass transport of surfactants and adsorption and desorption mechanisms at the interfaces need to be taken into account to include dynamic aspects.Extensive laboratory work and data extraction for different reactor geometries, as well as for more superficial velocities and DDA concentration levels, could provide a solid statistical knowledge base on these relationships.An empirical approach could be applied to the additional data with the aim of establishing correction factors for the numerically estimated k L a values for pure water.

Conclusion
We have highlighted the capabilities of the EE mixture model in terms of suitably capturing the physics of the aeration process.More specifically, we have presented the limitations of using Higbie's penetration theory at higher DDA concentrations.Research on models for k L and resistance to bubble movement for aqueous surfactant solutions appears to be somewhat limited.Moreover, to model the flow and oxygenation in surfactant solutions, further studies are recommended with the aim of developing theoretical models with support from experimental measurements to describe how these parameters are affected by the surfactant concentration.Due to the limitations of how the EE mixture model handles bubble dynamics at higher DDA concentrations, an empirical approach could be utilized by adding extensive laboratory data on the aeration process with the aim of establishing correction factors for the numerically estimated k L a values for pure water.The main conclusions drawn from the results of the study can be summarized as follows: i.The model can reliably predict the oxygen-transport rate in a bubble column with pure water at superficial velocities of 0.24 and 0.47 cm/s, but for higher superficial velocities, the model is less suitable as the flow regime becomes more heterogeneous and a more advanced two-phase model may need to be applied, such as the EE model, whereby a momentum equation for both faces is considered instead of one single face for the mixture.ii.To fully validate the model, the underlying theory, and the methodology, it is recommended that the model is tested to simulate systems with other reactor and air disperser geometries.iii.The entire model may be further developed according to one of two identified perspectives: an empirical perspective, where correction factors for primarily the mass-transport coefficient are developed and applied, and a theoretical perspective, where other mathematical models are used that are better suited for surfactant solutions.To do so, extensive laboratory work and data extraction on different reactor geometries should be conducted, along with more measurements at different superficial velocities and DDA concentrations in order to establish a solid statistical knowledge base on these relationships.iv.The flow resistance relationships in the EE mixture model are somewhat limited for non-uniform bubble dynamics.
Euler-Lagrange models (level-set, phase-field, or volume of fluid), which emphasize detailed analysis of single bubbles, could provide a better understanding of how to model flow resistance dynamics, i.e., using multiscale modeling principles to improve the current EE mixture and advection-diffusion models.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Boundary conditions for the computational domain, where the black text indicates the boundary conditions for the EE mixture model and the orange text represents boundary conditions for the advection-diffusion model.The dashed line cuts off the reactor and illustrates outlet conditions.

Fig. 2 .
Fig. 2. Example of a photograph of the reactor during operation, and dimensional measurement of the bubbles in AutoCAD (superficial velocity 0.71 cm/s, DDA concentration 9.375 mg/L).

Fig. 3 .
Fig. 3. Comparisons of DO concentrations in the reactor over time between experimental data ("exp") and COMSOL-simulated results for pure water ("mod").

Fig. 4 .
Fig. 4. Linear regressions of experimental and simulated data on the mass-transfer coefficient k L a as well as comparison of fitted curves of volumetric mass-transport coefficient k L a for experimental and simulated runs.

Fig. 5 .
Fig. 5. Linear regressions of experimental and simulation data on air volume fraction ϕ g as well as comparison of fitted curves of air volume fraction ϕ g for experimental and simulated runs.

Fig. 6 .
Fig. 6.Average slip velocities between air bubbles and water (u slip ) at different combinations of superficial velocities and DDA concentrations, the dashed lines are regression curves.

Fig. 7 .
Fig. 7. Comparison of the volumetric mass transfer coefficient, the gas hold up and the average slip velocity for the 2-dimensional axial-symmetric and fully 3-dimensional domain.

Fig. 8 .
Fig. 8. Contour plot of DO close to the inlet at the start of the simulation run as well as contour plots of air volume fraction, slip velocity, and volumetric mass-transfer coefficient close to the inlet at the start of the simulation run.

Table 1
Comparison of mass-transfer coefficients from experimental and simulated results as well as of air volume fractions from experimental and simulated results.