University of Birmingham Towards a robust CFD model for aeration tanks for sewage treatment – a lab-scale study

Detailed insight into the hydrodynamics of aeration tanks is of crucial importance for improvements intreatmentefficiency,optimizationoftheprocessdesignandenergy-efficientoperation.Thesefac-torshavetriggeredincreasinginterestintheuseofComputationalFluidDynamics(CFD)toevaluateperformanceofwastewatertreatmentsystems.Whilstfactorssuchasincorrectinputassumptions,poormodelchoiceandexcessivesimplificationshavebeenrecognizedaspotentialsourcesofout-puterrors,thereremainsaneedtoidentifythemostrobuststrategytofaithfullysimulateaerationtankperformance.Therefore,thefocusofthisworkwastoundertakerigoroustransientsimula-tions of the hydrodynamics and oxygen mass transfer in a lab-scale aeration tank in order to work towards the development of robust modeling guidelines for activated sludge systems. Unlike most previous CFD analyses of aeration systems, the work reported here employed the shear stress trans- port (SST) k − ω turbulence model to account for the turbulent interactions between the phases inducing bubble breakup and/or coalescence, and as a consequence, promoting the formation of bubbles of different sizes and shapes. The results obtained were compared with those arising from an analysis using the standard k − ε ( sk − ε ) model – and assuming fixed bubble diameter- the most common CFD modeling framework used within the wastewater modeling community. Model vali- dation was achieved using acoustic Doppler velocimetry and particle image velocimetry techniques, and experimentally derived oxygen mass transfer data. Limitations of both turbulence models used and modeling assumptions concerning bubbly flow are discussed. The benefits of the SST k − ω turbulence model are demonstrated, but the need to balance the increased computational expense of this approach compared to the sk − ε model and, indeed, bubble flow modeling are recognized. Thus, this paper presents the first rigorous analysis of turbulence model and bubble flow generation models together for activated sludge system optimization. Abbreviations: ADV: Acoustic Doppler velocimetry; AS: Activated sludge; ASM1: Activated Sludge Model No. 1; BOD: Biochemical oxygen CCD: Charge-coupled device; Computational fluiddynamics;COD:Chemicaloxygendemand;CPU:Centralprocessingunit;DO:Dissolvedoxygen;GCI:Gridconvergenceindex;HPC:Highperformancecomputing;IAC:Interfacialareaconcentration;MRF:Multiplereferenceframe;MLSS:Mixedliquorsuspendedsolids;PBM:Populationbalancemod-els;PIV:Particleimagevelocimetry;PST:Phase-spacethresholding;RAM:Randomaccessmemory;RANS:Reynolds-averagedNavier-Stokes;SNR:Signal-to-noise-ratio;SST:Shearstresstransport Eulerian model;


Introduction
Aeration is an essential element of the activated sludge (AS) process in many biological wastewater treatment plants, but it is also the most energy intensive, accounting for up to 75% of the total energy expenditure (Reardon, 1995;Rieger, Alex, Gujer, & Siegrist, 2006). AS systems require enhanced transfer of oxygen to maintain the biodegradation and nitrification processes, and hence the dissolved oxygen (DO) concentration in aeration tanks is a key process variable controlling biological oxygen demand (BOD) and nutrient removal efficiencies, as well as the overall operating cost of the process. Therefore, when considering the complex behavior of typical AS systems (Karpinska & Bridgeman, 2016;Karpinska, Dias, Boaventura, & Santos, 2015;Pereira et al., 2012), a reliable forecast of oxygen transfer, based on detailed insight into the hydrodynamics of the aeration tanks, is crucial for optimization of the process design and operation and for improvement of the biochemical conversion efficiencies. The recent increase in computational power and commercial availability of advanced software packages intended for the solution and visualization of complex flow problems found in industrial processes has contributed to the successful spread of computational fluid dynamics (CFD) within both academia and industry (Kochevsky, 2004). In particular, over the last two decades, continuous growth of the application of CFD to wastewater processes has been observed (Karpinska & Bridgeman, 2016). However, complete CFD simulation of multiphase flow in an AS tank coupled with biokinetic reactions remains a challenge due to the complexity of the models involved and the solution accuracy demanding a high level of mesh refinement, resulting in major random access memory (RAM) and central processing unit (CPU) requirements and long computational run times (Karpinska et al., 2015). Consequently, common practice is to simplify the modeling approach and simulate individual components of the AS system separately in a computationally efficient manner, and to couple the results afterwards (Pereira et al., 2012). Consequently, the majority of CFD studies concerning AS tank performance have focused on only one aspect, and have usually employed the most computationally inexpensive modeling scenario, based on steady-state or, occasionally, a transient approach and the sk − ε turbulence model. However, oversimplification of the CFD model and a lack of calibration will lead to errors in the simulation results and validation data. Moreover, despite the vast opportunities for the use of CFD modeling in wastewater engineering, a lack of educational support, training and targeted guidelines has been recognized as a core reason for frequent software misuse by inexperienced users (Wicklein et al., 2016), usually originating from incorrect problem statement, poor model choice, incorrect setup assumptions or misinterpretation of the results (Nopens et al., 2012).
Unambiguous guidelines for CFD modeling of aeration tank have not been defined so far. Thus, there is a need to elucidate a robust, 'core' multiphase model that will yield reliable results in terms of accurate local oxygen transfer rates -this being a critical parameter for prediction of the biokinetic reaction rates. Consequently, the objective of the work presented in this paper is to extend previous AS modeling work and to develop a transient CFD model of a lab-scale aeration tank, where the gas-liquid flow is modeled using a Eulerian model and the turbulence is simulated with two different dispersed models. Model outputs are validated against experimental results. Modeling scenarios initially assume fixed bubble size. Subsequently, the modeling scheme is enhanced via the introduction of bubble interaction models accounting for bubble coalescence and breakup. The impact of the turbulence and bubbly flow models and the imposed bubble size on the flow field, induced by the aeration and mixing systems, gas holdup and mass transfer, are also evaluated. The outcomes are compared with the results of a simpler approach (which is commonly used among wastewater modelers) and with experimental data, facilitating identification of the most accurate model, which is able to reproduce hydrodynamics and mass transfer in the aeration tank.

Previous modeling approaches
Modeling of aeration tanks based on a neutral density approach with steady Reynolds-averaged Navier-Stokes (RANS) simulations and the sk − ε turbulence model enables a relatively fast and straightforward prediction of the liquid-gas velocity fields and gas holdup; i.e. the percentage by volume of the gas (air) in the multiphase mixture in the reactor. Thus, this approach has been the preferred method by which to evaluate the performance of aeration and mixing systems aimed at process optimization through 'tune-to-benefit' operating parameters, particularly as it enables faster setup of lab-scale validation, involving the use of tap water and air only. However, as the suspended phase is neglected in this approach, the use of neutral-density models may lead to over-prediction of the degree of mixing  when translating the CFD results directly to fullscale agitated AS systems operating at high mixed liquor suspended solids (MLSS) concentrations, and where the occurrence of regions characterized by low liquid velocities has been identified. Therefore, work has been undertaken on the development of numerical tools based on a Eulerian twofluid model, which is able to assess the impact of the gas-liquid hydrodynamics on the global and local mass transfer coefficients in pilot-and full-scale oxidation ditches (one of several well-known modifications to the AS process) aerated with diffusers and agitated with slowspeed mixers (Cockx, Do-Quang, Audic, Liné, & Roustan, 2001;Do-Quang, Cockx, Liné, & Roustan, 1998;Fayolle, Cockx, Gillot, Roustan, & Héduit, 2007). The results from the simulations based on the fixed bubble size show that the gas holdup and oxygen transfer in closed-loop AS tanks depend on the axial velocity of the fluid imparted by the mixers. Moreover, accurate prediction of the overall mass transfer coefficient, k L a, relies on correct estimation of the initial bubble size at the diffuser level, requiring either in situ bubble diameter measurements (Fayolle et al., 2010) or implementation of an add-on model to estimate inlet bubble size. Similarly, the need for bubble diameter calibration has been emphasized in more recent CFD work by Terashima, So, Goel, and Yasui (2016), involving use of the same modeling approach to determine k L a values in a number of full-scale activated sludge systems and clean water tanks that differed in their dimensions, diffuser types (coarse and fine-pore, ceramic, plastic and membrane diffusers), their configuration (single and dual spiral roll) and operating airflow rates. Others (Cockx et al., 2001;Talvy, Cockx, & Line, 2007) have explored the use of a Eulerian model to determine axial dispersion and oxygen mass transfer in a pilot-scale airlift reactor. The simplified simulation strategy used assumed that the multiphase flow in the reactor was hydrodynamically in steady state, and thus only the transport equation for a source term representing oxygen recovery from air bubbles was solved. While numerical results obtained by Cockx et al. (2001) were relatively close to the experimental data, the discrepancies in values of k L a observed in the work by Talvy et al. (2007) were explained by the variable concentration of oxygen in bubbles assumed to be of constant size in simulations, and nonuniform mass transfer in the reactor.
In more recent work, Yang et al. (2011) developed CFD simulations based on a mixture gas-liquid model, which were used to optimize oxygen transfer in a fullscale multichannel oxidation ditch equipped with surface rotors and submerged mixers. Additional transport equations representing oxygen transfer and BOD removal were implemented. The results from this work showed that while computationally inexpensive RANS simulations can be used for design studies, the validity of this approach, based on averaged flow and concentration fields for optimization of the operating parameters, is rather limited. Gresch, Armbruster, Braun, and Gujer (2011) analyzed the flow patterns induced by porous diffusers in a rectangular, activated sludge tank. Similar to previous work, the fluid flow was simulated with a Eulerian twofluid model, and a fixed bubble diameter was defined for the dispersed phase. Unlike previous studies, the physical properties of the continuous phase were approximated to those of activated sludge, while the SST k − ω turbulence model was employed. Although oxygen transfer was not modeled, an additional transport equation for ammonia removal was solved. The outcomes from the CFD analysis showed that in 'mixerless' plug-flow aeration tanks, the flow field (as a function of diffuser layout) plays a significant role in air holdup and intensity of axial mixing.
A small number of researchers have implemented full, three-phase flow coupled with biological reactions, as described by the Activated Sludge Model No. 1 (ASM1) (Henze, Gujer, Mino, & van Loosdrecht, 2000) to simulate different AS tanks. Le Moullec, Gentric, Potier, and Leclerc (2010) studied a lab-scale AS tank aerated with a porous tube using a Eulerian two-fluid model enhanced with species transport, with source terms defining oxygen transfer, biokinetic reactions and biomass content. The simulations were based on the steady-state RANS approach and an assumption of fully soluble floc phase and the fixed bubble size. The chemical oxygen demand (COD) concentration profile predicted by the CFD model was in good agreement with the experimental values. However, overestimation of the simulated global k L a value and nitrate concentrations, and underestimation of the ammonia content, were observed and associated with excessive model simplifications necessitated by the feasibility of the simulation runs in realistic timeframes. More recently, an integrated hydrodynamic-biokinetic model calibrated experimentally was implemented to study a full-scale, closed-loop Carrousel-type bioreactor (Rehman, Maere, Vesvikar, Amerlinck, & Nopens, 2014). Contrary to the earlier CFD-ASM approach, the simulations considered the use of an algebraic-slip mixture model and the realizable k − ε turbulence model, and accounted for the bulk density computed from the suspended solids concentration. The resulting velocity field and concentration maps of DO and ammonium concentration allowed quantification of inhomogeneous mixing, and provided useful data in the context of reaction rates, energy efficiency and process control. Similarly, Lei and Ni (2014) used a calibrated CFD-ASM approach, based on a mixture and sk − ε model to represent fluid motion, phase interactions and biokinetics in a pilotscale Carrousel ditch equipped with mechanical aerators. The CFD model included solids transport and settling assuming pseudo-solid properties of the sludge flocs. The model was shown to be able to represent correctly not only the trends in hydrodynamics, oxygen transfer, and biological pollutant removal, but also the kinematics of sludge settling.
Further, comprehensive discussion on the CFD models applied to analyze, evaluate and optimize performance of the existing AS systems can be found in Karpinska and Bridgeman (2016) and Samstag et al. (2016).

Lab-scale aeration tank
Experimental studies were performed in a laboratory scale, two-phase, liquid-gas reactor, 0.50 × 0.10 × 0.26 m (length × width × depth) with designed active volume of 10 L (Figure 1). The tank was fitted with a fine-pore aeration system, comprising two diffuser tubes extending along the tank bottom, supplied with atmospheric air. This arrangement was designed to simulate a full floor grid configuration (i.e. two rows of diffusers) rather than a dual spiral roll layout. A peristaltic pump with a digital rpm speed controller (Watson Marlow Sci-Q 300) ensured continuous operation of the aeration tank at constant liquid flow rate of 0.03 m 3 d −1 . The aeration system was operated at varying flow rates, ranging from 0.1 to 1.2 L min −1 , regulated by means of two rotameters (Platon GTF1AHD). The tank content was homogenized by means of two overhead vertical shaft impellers (Heidolph RZR 2021) placed above the diffusers with straight blades of dimensions 0.08 × 0.04 m and 0.08 × 0.05 m. The operating scenarios considered are outlined in Table 1.

Velocity field measurement
The choice of an experimental method for robust determination of the velocity field in the lab-scale aeration tank was complicated by the mutual interference of the phases during the measurements. An acoustic Doppler velocimetry (ADV) probe was used to determine the 3D velocity of the liquid phase, whilst the nonintrusive laser particle image velocimetry (PIV) method was considered to be the most suitable method to assess air phase velocity field.
Acoustic Doppler velocimetry. The operation of an ADV probe placed in water is based on the emission of the short acoustic pulse at high frequency by a transmitter beam. The pulse traveling through the water column is scattered by the particles suspended in the water and reflected back to the ADV's receiver. The frequency of the Doppler shift between the transmitted and reflected acoustic pulses allows for calculation of the instantaneous fluid velocities in the measurement point. While the use  of ADV is straightforward, the measurements conducted in a bubbly flow environment, such as in aeration tanks, are difficult to analyze. The recorded raw ADV velocity data contain noise related mainly to the presence of the dispersed air bubbles, which move with different velocities compared to the liquid phase, and, to a lesser extent, due to Doppler signal aliasing (Mori, Suzuki, & Kakuno, 2007). Consequently, the output ADV velocity timeseries containing spikes of invalid data require despiking based on filtering algorithms and spike replacement methods (Jesson, Sterling, & Bridgeman, 2013). In this work, the velocity of the liquid phase was measured using Nortek Vectrino Plus ADV, operating at the acoustic frequency of 10 MHz and sampling rate of 200 Hz. The vertical extent of the sampling volume was 7 mm. The ADV probe was submerged at 4 cm below the fluid surface and aligned with the central cross-section along the tank. The velocity was measured in three locations, tagged as P1, P2 and P3, being the middle cross-section through the tank (between the impellers' shafts) and 6 cm from the front and back walls. For each operating condition listed in Table 1, a total of 120,000 velocity data points were collected. The raw ADV data were processed using the Velocity Signal Analyser software developed at the University of Birmingham (Jesson et al., 2013;Jesson, Bridgeman, & Sterling, 2015). The velocity time-series were despiked using correlation and signal-to-noiseratio (SNR) pre-filter, modified phase-space thresholding (PST) despiking filter (Parsheh, Sotiropoulos, & Porté-Agel, 2010) and linear interpolation method of spike replacement.
The first method is based on specifying an acceptable limiting value of signal correlation and SNR (recorded by the instruments), below which the value is identified as a spike and excluded from the analyzed data set (Jesson et al., 2015). The lower limits for correlation and SNR of 70% and 20 dB were set, according to manufacturer recommendations. In the modified PST filtering method, the raw data series are subject to a preconditioning step, which defines the upper limits for the values, allowing detection of the large magnitude spikes, and the lower limit for the valid data points in the vicinity of the large spikes, which should remain unchanged after despiking. This procedure is followed by application of the standard PST filter (Goring & Nikora, 2002) -an iterative technique using the minimal volume closing ellipsoid method to determine the statistical limit of valid fluctuating velocity components and their first-and secondorder time derivatives. The values falling outside the ellipse are labeled as spikes and removed from the data set. The data points removed by the despiking procedure are reconstructed using linear interpolation between the valid values on either side of the spike (Jesson et al., 2015).
Particle Image Velocimetry. The laser source used in the experiments was a Nd:YAG double-head pulsed laser (Litron Lasers, Model Nano L 50-100 PIV), emitting 4 ns pulsed beams with maximum energy of 400 mJ/pulse at a maximum rate of 100 Hz and a wavelength of 532 nm. The laser sheet was transported by means of an articulated optical arm with internal mirrors (TSI LaserPulse TM Light Arm, Model 610015) and the light sheet optics were positioned to illuminate the flow in the vertical cross-section throughout the tank along the diffuser's axis.
For each operating condition in the aeration tank, laser-synchronized, frame-straddled images of the illuminated bubbly flow were captured with a TSI PowerView TM Plus 4MP CCD camera, with resolution of 2048 × 2048 pixels and 12-bit output. A synchronizer (TSI LaserPulse TM , Model 610035) was used to control the timing between laser pulses, charge-coupled device (CCD) camera and image acquisition.
Each image captured a bubbly flow region of 0.25 × 0.20 m, corresponding to one half of the illuminated plane through the tank. Hence, in order to assess the average air flow patterns in the entire plane through the reactor, measurements were made for three camera positions aligned with the impellers' shafts and the center of the tank, thus capturing 'left' and 'right' near-wall regions, impellers and the central section of the tank. Spatial resolution of the measurements was 98.5 μm pixel −1 . The image frames captured by the CCD camera were transferred to a PC via a frame grabber.
TSI Insight 4G TM software was used for set-up, image acquisition and processing of the batch of raw images. For each experiment, 100 image pairs were recorded. Data sets composed of instantaneous vector fields in the selected regions-of-interest obtained with Insight 4G were further processed using Tecplot Focus 2013 software, yielding average velocity vector maps.

Oxygen transfer measurements
Measurements of the oxygen transfer in clean water and for the different operating conditions listed in Table 1 were carried out using a standard unsteadystate, reaeration method (ASCE, 2007;Von Sperling, 2007). Anhydrous sodium sulfite, Na 2 SO 3 , and crystalline cobalt (II) chloride hexahydrate, CoCl 2 ·6H 2 O, were used as deoxygenation chemicals. The experimental procedures and techniques used in preparation of the solutions, deoxygenation of the test water and hydraulic stabilization of the aeration system were performed according to the ASCE Standard (ASCE, 2007).
DO concentration and water temperature were continuously controlled with a portable dual channel multimeter (Hach Lange HQ40D) equipped with digital luminescent DO probes (IntelliCAL TM LDO10103), placed in the tank and in a transparent flow-cell adapter inserted into the outflow tubing from the tank. In view of the aeration tank size, distribution of mixers and aeration performance, the outflow was assumed to be representative of the total reactor content. During the first few minutes of the reaeration process, values of DO concentration and temperature in the measurement point were recorded every 0.5 s. Subsequently, the time interval for data acquisition was gradually increased by a few seconds, up to ca. 2 min. To minimize the occurrence of discrepancies in the results, the reaeration test was conducted three times for each operating scenario and for experimental conditions as close as possible to the standard (zero salinity, water temperature 20°C and pressure of 1 atm).
k L a values were determined using the log-deficit method (Mueller, Boyle, & Pöpel, 2002;Stenstrom, Leu, & Jiang, 2006), expressed as: where k L a is volumetric mass transfer coefficient, C * L denotes oxygen saturation concentration, C 0 is initial DO concentration at time t = t 0 , C L is the concentration in the bulk liquid phase, and the terms C * L − C 0 and C * L − C L represent the degrees of under-saturation at time t 0 , and after time t, respectively.
The value of k L a determined from the reaeration tests was corrected to standard conditions (i.e. 20°C).
where k L a T denotes k L a coefficient at any temperature T, k L a 20 is the k L a coefficient at the standard temperature of 20°C, and θ is the temperature coefficient. The recommended value of θ is 1.024 (Stenstrom & Gilbert, 1981).

3D geometry
The 3D geometry of the aeration tank was designed using ANSYS 15.0 Design Modeler preprocessor. The dimensions and placement of the inlets and outlet correspond to the lab-scale tank layout. Diffusers were designed as parallel, narrow, solid structures (0.46 × 0.015 × 0.015 m) located 20 mm from the lateral walls with the impellers placed 25 mm above the diffusers.

Hydrodynamics and mass transfer
Simulations of the hydrodynamics and mass transfer in the lab-scale aeration tank were performed using ANSYS 15.0 Fluent CFD software. Each simulation was run in parallel on the University of Birmingham Blue-BEAR Linux High Performance Computing (HPC) Cluster using dual-processor 8-core (16 cores/node) 64-bit 2.2 GHz Intel Sandy Bridge E5-2660 worker node with 32 GB of RAM. Bubbly flow in the lab-scale aeration tank was simulated with a Eulerian two-fluid model derived from RANS equations. The governing equations representing conservation of mass and momentum for continuous and dispersed phases have been described in detail elsewhere (Azzopardi et al., 2011;Karpinska & Bridgeman, 2016;Ratkovich, 2010). Two different two-equation turbulence models were considered, specifically the sk − ε model (Launder & Spalding, 1974) with scalable wall functions, and the shear stress transport (SST) k − ω model (Menter, 1994). A comprehensive review of these models, their applicability and limitations can be found in Bridgeman, Jefferson, and Parsons (2009) and Karpinska and Bridgeman (2016).
The drag coefficient used in the work reported here was obtained using Grace et al.'s model (Clift, Grace, & Weber, 1978). This model significantly improves upon commonly used drag models (i.e. Schiller Naumann and Morsi Alexander, which are based on the assumption of rigid-sphere bubbles), by adjusting the drag coefficient over a wide range of shapes of bubbles (spherical, ellipsoidal and capped). The Brucato correlation (Brucato, Grisafi, & Montante, 1998) was used to modify the drag factor, accounting for the impact of the continuous phase turbulence. In order to account for the dispersion of the secondary phase due to the transport by turbulent fluid motion, an additional turbulent dispersion force was modeled using the Simonin approach (Simonin & Viollet, 1990). The drift velocity term, which originated from Tchen-theory correlations for the dispersed two-equation models (Hinze, 1975), was able to predict the interfacial turbulent momentum transfer. Moreover, the turbulent interactions between the dispersed and continuous phases were accounted for via implementation of additional source terms in transport equations for k and ε for the continuous phase, using the model proposed by Simonin and Viollet (1990). The governing equations related to the above models are included as Supplementary Information.
The transport equation for the interfacial area concentration (IAC) accounting for coalescence and breakage effects is defined as (Ishii & Kim, 2001): where a i is the IAC, v i is interfacial velocity, α G is volume fraction of secondary (gas) phase and v G is its velocity, S RC and S WE are the coalescence sink terms due to random collisions and wake entrainment, S TI represents brakeage source term due to turbulent impact and φ Ph accounts for the contributions to the interfacial area concentration due to the phase change as a result of nucleation, evaporation or condensation. In the case of adiabatic, isothermal two-phase flows, the term φ Ph is neglected.
The shape factor, ψ, is defined as: where d bS and d e are the Sauter-mean and volume equivalent diameters.
In the work reported here, the effects of bubble breakage and coalescence on the hydrodynamics in the aeration tank were modeled using the Hibiki-Ishii approach (Hibiki & Ishii, 2000). The IAC model assumes that the bubbles behave like ideal gas molecules (Coulaloglou & Tavlarides, 1977). The coalescence is caused by random bubble collisions induced by the turbulence in the liquid phase. The coalescence rate equivalent to the sink term S RC in Equation (3) is given by: where f C is collision frequency; n b is bubble number density; λ C is coalescence efficiency; C is empirically determined adjustable value, which for bubbly flow equals 0.188; ε L is turbulent kinetic energy dissipation rate of the continuous (liquid) phase; α max is maximum allowable air fraction; d b denotes bubble diameter; ρ L is the density of the primary phase and σ is surface tension between the phases; and K C is a model constant ( = 1.29).
The bubble breakup occurs as a result of the collision between the turbulent eddy and the bubble. The bubble breakup rate, expressed as S TI is obtained via: where f B is bubble-eddy random collision frequency, n e is number of eddies of wave number per volume of twophase mixture, λ B is breakup efficiency, B is empirically determined adjustable value, which for bubbly flow is 0.264, and K B is a model constant equal to 1.37. The interfacial oxygen mass transfer occurring between air bubbles and liquid is: where L L is interfacial mass transfer between air bubble and liquid, k L is the local mass transfer coefficient, a is the interfacial area and the gradient (C * L − C L ) is the driving force causing oxygen transfer.
The interfacial area, a, is calculated as: The local mass transfer coefficient, k L , is obtained from the Higbie penetration theory (Higbie, 1935): where D L is the molecular diffusion coefficient of oxygen in water at 20°C and μ L denotes the dynamic viscosity of phase L.

Boundary conditions and simulations scheme
For the purposes of the study recorded here and validation techniques requiring the use of clean water, modeling of hydrodynamics and mass transfer in the lab-scale aeration tank was based on gas-liquid neutral density simulations. Therefore, water with a constant density of 998.2 kg m −3 and dynamic viscosity of 0.001 Pa s was set as a continuous phase. Air having density of 1.225 kg m −3 was defined as a secondary phase.
The assumed average hydraulic retention time was 8 h (a typical value for biological nutrient removal process tanks, such as plug flow or step-feed AS systems [Jenkins & Wanner, 2014]). Velocity inlet boundary conditions were imposed on the circular surfaces representing influent ports, and the local velocity of the water in the normal direction to the boundary was set as 0.006 m s −1 . A pressure outlet boundary condition was set on the outflow, with backflow of the air phase fraction of 0. A velocity inlet boundary condition was imposed on the active surfaces of the diffusers, with local air velocity for each diffuser corresponding to the air flow rates of 0.1 × 10 −3 to 0.5 × 10 −3 m s −1 . A degassing boundary condition was imposed on the fluid surface.
The multiple reference frame (MRF) approach was used to model impellers operating at two speeds: 50 and 150 rpm, and a no-slip condition was imposed on the lateral walls and the bottom of the aeration tank. A movingwall boundary condition was set on the surfaces of the impeller blades and shafts.
During the simulations, the operating pressure set at the water surface was 101325 Pa and the flow was assumed to be isothermal (operating temperature of 293 K). Acceleration due to gravity was set at 9.81 m s −2 and the specified operating density was set to that of the dispersed phase -i.e. 1.225 kg m −3 .
Convergence criteria for the solutions were 10 −6 . For the sake of stability of convergence of the SST k − ω model, the first 10 5 iterations were run with the initial time step size ( t) of 0.0001 s. With stable residual monitors, t was increased to 0.01 s. After the first 7000 s of the flow time, negligible changes in simulation outputs were observed. Nevertheless, all the simulation runs captured a flow time of 14,400 s (4 h).
The simulation schemes used to study the lab-scale aeration tank are outlined in Table 2.

Mesh
The fluid domain was divided into groups of bodies with either sweepable or unsweepable topologies. The computational grid generated through the sweeping algorithm was based on the projection of the quad-predominated source surface mesh created for the top wall onto the target surface within the considered body or group of bodies. Two features of the source surface mesh were enabled, viz. inflation of the lateral edges of the rotating zone, adjacent to the impeller and shaft) and the sweep bias (Figure 2). The remaining unsweepable bodies were meshed using the patch conforming method, based on the triangular surface meshing algorithm. The face sizing function was enabled to allow mesh refinement in the inlet and outlet zones. The different features of the source mesh and predefined minimum element size  yielded four meshes composed of 1.5 to 3.3 million hexaand quadrilateral cells (Table 3). The Grid Convergence Index (GCI) approach (Roache, 1998a(Roache, , 1998b was employed to assess mesh density. The GCI is a recommended uncertainty estimator method (Celik & Karatekin, 1997) for uniform reporting of grid refinement studies outcomes through the prediction of a discretization error in the finest mesh relative to the converged numerical solution. Details of the GCI calculations and the results obtained are provided as Supplementary Information.
The results indicated that Mesh 3 was appropriate for subsequent modeling work.

Experimental validation of the CFD model
The main objective of the work presented in this paper is to identify the most reliable transient CFD model to simulate aeration tanks. As a result, the emphasis was on the model capability to reproduce experimental results obtained for a range of operating conditions. Figure 3 shows samples of raw ( Figure 3a) and clean ( Figure 3b) velocity time-series data in one direction (vertical component) collected during a 9-min period of bubble-influenced flow, from one measurement point in the lab-scale aeration tank. The clean data were obtained after 10 iterations of the despiking procedure.

ADV
The filtered ADV velocity time-series measured at three locations within the tank (P1, P2 and P3) in varying operating conditions are presented as Supplementary Information ( Figure SI 2-4). The contour maps of the water velocity in the cross-section through the middle of the tank obtained from the CFD simulations with different turbulence models and for a constant bubble size are presented in Figure 4.
For the constant mixing speed and airflow rates from 0.1 to 0.8 L min −1 , the average values of the water velocity measured with the ADV (Figures SI 2-4 a-c) are slightly above zero. These results are in close agreement with the outcomes of both turbulence models in   the zones corresponding to the points P1-P3, as seen in Figure 4a-c. The increase of the mean water velocity was observed in ADV data recorded at points P1 and P2 for the airflow rate of 1.2 L min −1 ( Figure SI 2-3  d). These changes were reproduced by the SST k − ω model, as seen in the evolution of the velocity contours and direction of the vectors (Figure 4d), but were clearly underestimated by the sk − ε model.
Despite visible changes in the vector and contour maps obtained with both CFD models, increases in water velocities recorded at locations P1-P3 for the 150 rpm mixing speed ( Figures SI 2-4 e) were underestimated by the sk − ε model (Figure 4e). However, these experimental results were reproduced well in the contour plots obtained with the SST k − ω model.
The ADV data were used as a benchmark for comparison with the outcomes of the CFD simulations. Figure 5 shows a bar chart representing differences between the measured and predicted velocity magnitudes by both models' water velocities at points P1-P3. While both turbulence models failed to reproduce the exact experimental values, the differences between the benchmark and CFD data obtained with the SST k − ω model for the higher airflow rates (0.8 and 1.2 L min −1 ) and the mixing speed of 150 rpm were distinctly smaller ( Figure 5). On the other hand, the sk − ε model performed well for the lower range of airflow rate (0.1-0.4 L min −1 ) and the mixing speed of 50 rpm, as seen for points P1 and P2 in Figure 5.
However, direct translation of ADV data to CFD results is not straightforward. One reason for the difference between measurement and simulation results is that both the sk − ε and SST k − ω models are based on the Bussinesq isotropic eddy viscosity assumption, which leads, especially in case of the sk − ε model, to inaccurate prediction of the flows driven by anisotropy of the normal Reynolds stresses and secondary shear stresses, and flows characterized by large extra strains, e.g. swirling flows (Bridgeman et al., 2009).
Moreover, in common with the vast majority of CFD studies of lab-scale systems, the simulations assumed ideal flow conditions in the tank but did not account for the interference from the submerged ADV sensor itself, which, considering tank dimensions, could influence to some extent the evolution of the flow patterns. Therefore, further analysis of the CFD and experimental results for the air phase is necessary for robust selection of the optimal CFD model. Figure 6 presents the air velocity vector maps in a vertical cross-section through the tank aligned with the diffuser's axis obtained from the CFD simulations using two different turbulence models for varying operating conditions overlying the corresponding PIV results. The CFD simulations shown considered fixed air bubble size at the diffuser level and the shape of the bubbles being dependent on the drag and turbulent interactions with the continuous phase. The corresponding pair of PIV figures represents the section of the analyzed flow field above the diffuser.

PIV
Regardless of the turbulence model applied, good agreement between the CFD simulation results and the measured air flow field was observed for a constant mixing speed of 50 rpm and air flow rates from 0.1 to 1.2 L min −1 (Figures 6a-d).  Nonetheless, a small number of differences between the simulated and measured data can also be identified. At 150 rpm mixing, the CFD studies yielded distinctly higher local velocity values than the PIV results in the region of the rotating blades ( Figure 6e). This difference may be explained by the fact that the applicability of the standard 2D PIV technique to bubbly flows is limited to low operating gas volume fractions, mainly due to the scattering of the light sheet on the bubbles' surface, causing distortions and biasing of the images (Sommerfeld, 2004). The adverse effects of the background noise of the image induced by the presence of the bubbles outside the measurement plane or changing bubble dimensions due to breakup/coalescence occurring on the surface of the rotating blades have been also recognized (Honkanen, Saarenrinne, & Larjo, 2003).
The PIV maps obtained for the impeller region (righthand side of the analyzed plane) indicated the formation of low-air-velocity zones adjacent to the rotating shaft, which is slightly underestimated by the CFD simulations (Figures 6a-c, e). The occurrence of small low-airvelocity areas just above the diffuser was also neglected by the CFD models. These discrepancies between the simulated and measured outcomes may be explained by the fact that the CFD models assumed ideal aeration conditions where the porosity of the diffuser does not change in time, producing a uniform air plume consisting of bubbles with identical 'inlet' diameters. Despite undertaking particular efforts to approximate experimental conditions to the simulated scenario (e.g. use of new, clean diffusers; carrying out the experiments in deionized water to avoid fouling; and minimizing flow perturbances by reducing air supply tubing to one per diffuser), maintaining uniform bubble distribution above the diffuser remained a challenge, as can be seen in Figure SI 5.
The choice of the airflow rate of 0.8 L min −1 for subsequent measurements was based on the visual inspection of several raw PIV images obtained for each of the operating conditions. The resulting air velocity field in the central cross-section and the CFD results obtained with different models are shown in Figure 7. It can be clearly seen that both numerical models reproduced the airflow field robustly, taking into account local velocity The results obtained in the PIV experiments confirm the validity of both CFD models enabled to simulate hydrodynamics of the two-phase lab-scale aeration tank.

Turbulence model comparison
The main objective of this section of the work was to select the most reliable transient CFD model able to reproduce experimental results obtained for a range of operating conditions, rather than optimization of the reactor performance. Accordingly, in this section, the impact of the turbulence model on the flow field and the gas holdup for varying operating conditions is discussed.

Impact of turbulence model on liquid phase velocity field
In two-phase aerated bioreactors, the overall mixing phenomena are linked to the dynamic interactions occurring between the ascending air plume introduced by the diffusers and the horizontal velocity of the liquid phase imparted by the momentum sources -that is, the impellers.
A contour map representing distribution of the water velocity vectors in the vertical cross-section through the middle of the aeration tank is shown in Figure 4. Regardless of the turbulence model used, a rotational speed of 50 rpm was insufficient to ensure effective mixing of the tank contents. In conditions of low airflow rates up to 0.8 L min −1 (Figures 4a-c), the occurrence of stagnant fluid regions, characterized by low or no mixing, is evident. A slight improvement in mixing was observed for the airflow rate to 1.2 L min −1 (Figure 4d). Here, the subtle but dynamic changes in vector field were captured particularly well by the SST k − ω model. Similar conclusions can be drawn when comparing velocity distribution in the horizontal cross-section through the tank at the level of the impellers' blades ( Figure 8). Increasing the airflow rate led to gradual disappearance of the stagnant zone between the impellers, as shown in contour maps obtained with the SST k − ω turbulence model (Figures 8a-d), but this was less pronounced in the case of the sk − ε model (Figures 8a-d). The effective use of the whole volume of the tank, due to increase of the water velocities within the analyzed cross-sections, is possible only for the mixing speed of 150 rpm (Figures 8e  and 4 e).
Differences between the results obtained with the two turbulence models are clear when comparing the water velocity patterns in the near-wall region and above the tank bottom (Figures 4 and 8). These differences are associated with the model sensitivity to solve the low-Re number boundary layer region at no-slip walls. It should be noted that although the viscous wall-bounded flow was not thoroughly studied in this work, the nondimensional distance from the wall to the first mesh node, y + , was in the range of 0.5 < y + < 10, allowing solution of the viscous sub-layer (y + < 5) and part of the adjacent buffer zone (5 < y + < 30) (Pope, 2000) using the SST k − ω model. However, the boundary layer region was not necessarily solved with sk − ε model, but was predicted via application of the appropriate wall functions, providing a numerical bridge between the viscosity-affected nearwall region and the inertia-dominated turbulent flow core.
Further differences between the models were observed in the evolution of water velocity patterns and vector fields in the region affected by rotation (Figures 4 and 8): these were dynamic and irregular for the SST k − ω model and more steady for the sk − ε model. These differences are a result of the poor performance of the sk − ε model in more complex flow scenarios characterized by, for instance, strong streamline curvature, vortices, rotating flows and the model's limited sensitivity to body forces due to rotation of the reference frame (Versteeg & Malalasekera, 1995). Further inaccuracies of the sk − ε model may originate from the entirely empirical ε equation (Pope, 2000). Another known drawback of the sk − ε model is the assumption of locally isotropic turbulence, where the prediction of the turbulent viscosity (μ t ) is based on the constant C μ ( = 0.09), whereas the SST k − ω enables modified μ t formulation with variable C μ to account for transport effects of principal turbulent shear stresses (Karpinska & Bridgeman, 2016). Figures 6 and 9 show the contour and vector maps representing air velocity in two parallel vertical cross-sections, which are extended along the diffuser and through the middle of the tank. For the cases shown, the CFD simulations considered fixed bubble size, with its shape dependent on the drag and turbulent interactions with the continuous phase.

Impact of the turbulence model on gas phase velocity field
The CFD results obtained with different models (Figure 6) demonstrated similar characteristics in distribution and direction of the velocity vectors in the analyzed plane aligned with the diffuser. Regardless of the turbulence model used, for the lowest airflow rate of 0.1 L min −1 (Figures 6a, 9a), the maximum values of the air velocities occurred in the region of the rotating blades.
While further increases of the operating airflow rate did not increase the local air velocities on the blades (Figures 6b-d and 9b-d), the formation of an air plume with higher velocities, especially in the region between impellers, can be discerned. Increasing the mixing intensity (Figures 6e and 9e) significantly increased the air velocity, most notably in the vicinity of the rotating blades and adjacent to the lateral walls.
Both models provided different predictions of the air velocities in the near-wall region and between the rotating impellers. The dynamic changes in the shape of the air plume between the impellers, and increase of the air velocities close to lateral wall, were reproduced by the SST k − ω model, but clearly underestimated by the sk − ε model (Figures 6a-e and 9a-e). These results support the findings from analysis of the water velocity field shown in Figures 4 and 8, concerning limitations of the sk − ε model and higher accuracy of the SST k − ω model in predicting the actual flow features in the viscous near-wall region and inertia-dominated rotating zones. It should be noted that the motion of the bubbles in the tank depends to a great extent on correct prediction of the liquid phase velocity and the turbulence quantities, k and ε, as they are enabled in the modeling of the interfacial turbulent momentum transfer, dispersion of the bubbles

Impact of turbulence model on gas holdup
Air holdup is one of the most important hydrodynamic parameters governing oxygen mass transfer efficiencies in aeration tanks. Figure 10 shows the air volume fractions in the vertical cross-section along the diffuser for varying operating conditions obtained with different turbulence models. It is known that air holdup depends on the superficial velocity of the air phase and the bubble diameter, and thus on the retention of air bubbles in the tank. Since all the simulation schemes considered the air phase as having a fixed bubble diameter of 0.001 m, it was expected that the air holdup would yield higher values with increasing operating air flow rates. Accordingly, regardless of the turbulence model applied, and under the same mixing conditions, the lowest airflow rate yielded the worst operating scenario and the lowest gas holdup of around 0.1% of the tank volume (Figure 10a). Increasing the airflow rate resulted in a proportional increase of the volumetric gas fraction (Figures 10b-d), with the maximum local values of the air holdup of 1.0% in the crosssection, distributed just above the diffuser. Increasing the mixing speed did not improve tank performance in terms of air holdup, as shown in Figure 10e. For the analyzed cross-section throughout the tank, the decrease in the air holdup is clearly pronounced in the rotating zones. This phenomenon is attributed to the increase in the superficial liquid velocity induced by the impellers (see Figure 4e), which neutralizes the vertical flow induced by the diffusers, resulting in reduced residence times of the air bubbles and their faster disengagement. Taking into account that the air volume fraction in the tank depends on the local air velocities (Figure 9), the results of the SST k − ω and sk − ε turbulence models differ in the distribution of local gas holdup in the mixing zones and near the wall (Figure 10). This is consistent with earlier discussions on the model's capacity to reproduce actual flow conditions in the agitated tank.

Bubbly flow models
The bubble diameter, the character and frequency of occurrence of interactions involving the dispersed and continuous phases are of crucial importance for the correct prediction of air velocities and holdup, and hence for assessment of the mass transfer phenomena occurring in the analyzed system. Figures 11 and 12 show contour maps of the air velocity and gas holdup obtained with the SST k − ω and sk − ε turbulence models for varying assumptions regarding bubbly flow -i.e. constant and varying bubble diameters from 0.5 to 1.0 mm (IAC model). The range of the bubble diameter was selected on the basis of a visual inspection of several picture frames captured by a CCD camera ( Figures SI 6 and 7). Two operating conditions are compared: air flow rate of 1.2 L min −1 and two mixing speeds of 50 rpm and 150 rpm. Considering a mixing speed of 50 rpm, the results obtained with the SST k − ω model accounting for the bubble interactions provided noticeably lower air velocities than did the simulation results obtained with constant bubble size, as seen in Figure 11a. However, the velocity and vector plots obtained in the cross-section aligned with the diffuser were in good agreement with the PIV data, reproducing the low-velocity region in the vicinity of the impellers and above the diffusers (Figure 13a). The occurrence of the lowest air velocities above the bottom (Figure 11a) indicates the presence of small bubbles with diameters of less than 1.0 mm, and consequently with lower-rise velocities. The higher water velocities on the rotating blades may contribute to the local breakup of the bubbles and their faster disengagement. Under these circumstances, lower air holdup compared to the case of a constant bubble size was expected, as shown in Figure 12a.

Impact of the bubble size on air velocity and gas holdup
The simulation of bubbly flow with the SST k − ω and IAC models for the ISO rpm mixing speed to 150 rpm yielded similar results to the constant bubble modeling scenario (Figure 11b), except for the local low air velocities above the rotating blades and the tank bottomwhich was associated with the release from the diffusers a swarm of bubbles with diameters smaller than 1.0 mm. The vector map of the air velocity in the cross-section through the diffuser was also in close agreement with the PIV data shown in Figure 13b. The effects of the bubble interactions are clear when comparing the contour plots of the air volume fraction in cross-section through the tank (Figure 12b). The air holdup values predicted by the IAC model were distinctly lower than the outcomes of the modeling scheme based on the constant bubble size. This outcome is likely to be a result of shortening the residence times of the bubbles in the tank. Increased mixing speed results in increased superficial water velocity acting on the plume of shredded bubbles with low superficial velocities, facilitating their dispersion and later release through the fluid surface. On the other hand, increased turbulence of the water near the rotating blades promotes random bubble collisions and formation of the larger bubble structures. The coalescing bubbles are characterized by higher-rise velocities (Figure 11b), resulting in their shorter residence time in the tank, and thereby in the low air holdup seen in Figure 12b.
Contrary to the results of the SST k − ω model, implementation of the IAC model did not yield significant differences in the air velocity and holdup maps obtained with the sk − ε model, as seen in Figures 11c-d  and 12c-d. Irrespective of the mixing speed, the model accounting for bubble interactions and varying bubble diameter resulted in the occurrence of the narrow region of lower air velocities from the bottom to slightly above the diffusers (Figures 11c-d), linked to the presence of the small bubbles. While it is unclear whether bubble breakup took place, formation of the low velocity region just below the impellers was most likely due to the bubble coalescence, promoting formation of the bubbles with uniform diameters of 1.0 mm (and hence striking resemblances in the velocity contours to the results obtained  with fixed bubble diameter). Consequently, both bubbly flow models yielded remarkably similar contour maps of the air holdup in the analyzed cross-section, as seen in Figures 12c-d. To exclude the impact of the narrow bubble size range (0.5-1.0 mm) on the outcomes of the sk − ε model, an additional simulation run for varying bubble sizes, from 0.5 to 3.0 mm, was performed. The resulting maps of air velocity in the cross-section through the tank combined with the outcomes of the previous IAC simulations are show in Figure 14a-b. Regardless of the upper bubble size limit, both simulations resulted in the formation of a similar size and shape low-air-velocity region extended above the tank bottom. Increasing the maximum bubble size limit caused a sudden and significant increase of the air velocities in the tank, as shown in Figure 14a. These results confirm the findings from simulations obtained for the bubble size of 0.5-1.0 mm related to the predominant role of coalescence in the bubble interactions being modeled. As coalescence leads to an increase of bubble size, shortening retention times; the resulting bubble holdup was lower than for the narrow-diameter size range (Figure 14c-d), but not as low as for the SST k − ω model (Figure 12a-b).
The outcomes obtained with the IAC and both turbulence models revealed substantial differences in terms of predicting the bubble interactions and resulting air holdups. The root cause is model accuracy in prediction of the local velocities of the continuous phase and the turbulence parameters, k and ε. The turbulent kinetic energy dissipation rate, ε, and bubble diameter, are used in the prediction of bubble coalescence or breakage rates (see Equations (5) and (6)), which play a decisive role in assessing the local air velocities, and thus air holdup.
To confirm the above conclusions, the values of the turbulence kinetic energy dissipation rate of the liquid phase obtained with different models are shown in Figure 15. It can be seen that the empirical values of ε enabled by the sk − ε model are approximately twice Figure 15. Values of the turbulent kinetic energy dissipation rate obtained with different modeling schemes. those for the SST k − ω model, which explains the differences in the extent of predicted bubble interactions.

Mass transfer coefficient
ADV and PIV results did not lead to an unambiguous selection of the turbulence model. However, an analysis of the hydrodynamics in the aeration tank using different turbulence and bubble interaction models demonstrated substantial differences in predicted air velocity fields and air holdups. While the air velocity field obtained with IAC models may be in close agreement with PIV figures, an additional selection method was still required to assess the validity of the modeling approach, and the oxygen mass transfer coefficient was selected for this purpose.
While there are many hydraulic and operational parameters affecting the performance of diffused aeration systems, the oxygen mass transfer is controlled by two key factors: the local mass transfer coefficient and the total interphase area of the air bubbles. The local mass transfer coefficient, k L , is associated with the turbulence quantities of the liquid phase, namely ε L (Equation (9)), while the interphase area depends on the contact time between the rising air bubbles and the liquid, and hence on their diameter, velocities and, therefore, on the air holdup. The consequences of the use of different turbulence and bubbly flow models are more pronounced when one compares the hydrodynamic parameters and the standard volumetric mass transfer coefficient k L a 20 summarized in Table 4. To facilitate comparison of the CFD data with the measurements' outcomes, the values of the oxygen mass transfer coefficient determined in several experiments are also included in the table. Sample results of a selected reaeration test conducted in the lab-scale aeration tank are shown in Figure 16.
Irrespective of the mixing speed, when analyzing data obtained using the sk − ε model for the airflow rates of 0.8-1.2 L min −1 (Table 4), it is observed that the IAC model for bubble sizes from 0.5 to 1.0 mm yielded approximately the same results as the SST k − ω and sk − ε models did for fixed bubble size, for superficial air velocities, holdups and interfacial areas. The standard mass transfer coefficients were also overestimated compared with the experimental data. As concluded previously (Figures 11-12c-d), the reason for this originates from the overestimation of effects of coalescence uniformizing bubble sizes towards the assumed upper size limit; in this case, 1.0 mm. The same phenomenon was observed in the results from the sk − ε and IAC models, with bubble sizes of 0.5-3.0 mm (Figure 14). Larger bubble sizes resulted in their lower concentration in the tank, leading to a sharp decline in values of the specific interfacial area, and thus mass transfer coefficients (Table 4). The outcomes of the SST k − ω and IAC models (Table 4) differ markedly from the results obtained with fixed bubble size. In this case, due to simultaneous bubble breakup and coalescence, lower air velocities and holdups (shown in Figures 11a-b and 12a-b) result in a decrease of the specific interfacial areas by a factor of approximately 2 (Table 4). Consequently, the resulting standard mass transfer coefficients are distinctly lower than those obtained assuming fixed bubble size.
To determine the validity of the CFD results, the average experimental k L a 20 (Table 4) was used as a benchmark for comparison with the outcomes of different modeling schemes. A bar chart representing differences in experimental and predicted values of k L a 20 is shown in Figure 17. In the case of simulations performed for the mixing speed of 150 rpm (Figure 17a), the results of the SST k − ω and IAC models were seen to be the most accurate in prediction of the k L a 20 . In this case, the decisive factor in model performance was associated with accuracy in prediction of the rotating velocity field and the turbulence parameter, ε, governing the occurrence and intensity of the bubble interactions and a local mass transfer coefficient, k L . The same characteristics were observed in the results obtained for a mixing speed of 50 rpm and operating airflow rate of 0.8 L min −1 (Figure 17c). Slightly less pronounced superiority of the outcomes of SST k − ω and IAC model are reported for the air flow rate of 1.2 L min −1 (Figure 17b). While one of the reasons may be the experimental value of k L a 20 (average value from three experiments), the other possibility may be linked with the bubbly flow regime and the influence of the empirical adjustable model parameters in Equations (5) and (6) ( C = 0.188 and B = 0.264) on the actual coalescence and breakup rates. It should be noted that the latter model parameters were determined from experiments carried out on the vertical, rounded tubes (Hibiki & Ishii, 2000;Ishii, Kim, & Uhle, 2002;Wu, Kim, Ishii, & Beus, 1998).
Considering outcomes of both turbulence models with fixed bubble size and for low airflow rates of 0.1 and 0.4 L min −1 (Table 4), the experimental values of k L a 20 were well reproduced by the SST k − ω model, whereas they were clearly overestimated by the sk − ε model, as shown in Figure 17d. While the assumption of uniform bubble size would appear to fit the dispersed bubbly flow regime, which is characterized by small variation in the sizes and shapes of the bubbles, the differences in the model outcomes are more likely to be linked to the prediction of turbulence in the water phase. However, since the bubbly flow regime was not explicitly determined experimentally in this work (e.g. via the standard imaging technique), the balance between the coalescence and shredding of the bubbles in the mixing zones (not modeled in this case) cannot be excluded. Therefore, for the case of the lab-scale aeration tank simulated with SST k − ω, the solution accuracy will depend on determination of the relationship between the lowest liquid and gas superficial velocities, for which implementation of the coalescence and breakup models is necessary.
While the values obtained with the SST k − ω model and models accounting for the bubble interactions offer improved prediction of the actual hydrodynamics and mass transfer in the tank, it is also clear that successful simulation of the coalescence and breakage effects with the IAC models depends on correct estimation of the input data, namely the physical properties of the dispersed phase. When taking into account the turbulent flow regime occurring in the aeration tanks, the frequency and intensity of collisions, breakups and deformations yield a large diversity in the shapes and sizes of the bubbles (Karpinska Portela, 2013;Shaikh & Al-Dahhan, 2007;Takács, 2005). Nonetheless, a tendency to simplify the complex modeling framework for aeration tanks through the assumption of nondistributed scalar properties of the phases, i.e. fixed bubble diameter, is still predominant (Karpinska & Bridgeman, 2016). An enhanced multiphase approach accounting for bubble size distributions within the aeration tank requires implementation of population balance models (PBM). A comprehensive assessment of several common solution methods for PBM, viz. the discrete class size method (Hounslow, Ryall, & Marshall, 1988), the standard method of moments (Randolf & Larson, 1971) and the quadrature method of moments (Marchisio, Vigil, & Fox, 2003) can be found in Bridgeman et al. (2009). It should be emphasized that although PBM-CFD coupling has been successfully exploited in the chemical engineering sector to study bubble columns, airlift and stirred bioreactors (Dhanasekharan, Sanyal, Jain, & Haidari, 2005;Morchain, Gabelle, & Cockx, 2014;Wang, 2011), its application in modeling of aeration systems with suspended solids is still uncommon and is limited to only a very few examples (Karpinska & Bridgeman, 2016;Nopens, Torfs, Ducoste, Vanrolleghem, & Gernaey, 2015).

Modeling issues -Computational cost
A common criterion for the selection of a modeling scheme is related to RAM and CPU usage and the computational expense. In the work reported here, the computational cost was increased via the refined mesh requirements of the SST k − ω model and the complexity of the modeling approaches involved in the simulations of the hydrodynamics and mass transfer in lab-scale aeration tank. The computational time required to simulate a flow time of four hours using the sk − ε model and using the simplest modeling scheme (Table 2) based on the MRF, two-fluid model and the constant bubble size approach was approximately five days using BlueBEAR HPC. For the SST k − ω model, which is known to be the most expensive of the RANS-based two-equation models, the same simulation flow time was achieved in eight days. These computational requirements in terms of RAM and CPU and computational time were approximately doubled for the simulation schemes accounting for bubble interaction (an extra transport equation for interfacial area concentration).

Lab-scale approach for full-scale AS tanks
Irrespective of the purpose of the numerical analysis, experimental calibration of the simulation input data and validation of CFD results is still required, especially in the complex flow situations found in full-scale AS tanks. Hence, the CFD analysis procedure performed for the lab-scale aeration tank presented in this paper can be used for troubleshooting full-scale aeration tanks. The procedure is based on downscaling the AS system being analyzed to the lab-or pilot-scale model, making it feasible to simulate with CFD codes using a Eulerian multiphase model, followed by validation using, for example, stereo-PIV based techniques. While use of flow-intrusive measurement techniques is limited due to the lab-scale system dimensions, a supporting selective validation method is highly recommended. For determination of the relevant flow characteristics, such as formation and shape of the bubble plumes in near-wall/bottom regions or analysis of the flow patterns induced by the mixer and resulting mass transfer coefficient, application of the computationally expensive SST k − ω model enhanced with models accounting for bubble breakage, coalescence and bubble size distribution is more appropriate. Furthermore, the proposed 'core' hydrodynamic model developed for water and air can be expanded through implementation of the calibrated submodel, adjusting density of the continuous phase to account for the presence of MLSS, thus extending its applicability for AS systems. The outcomes are expected to provide a robust prediction of the DO concentrations in the tank, which is of crucial importance for treatment performance evaluation (COD and ammonia removal) and for optimization studies based on the controller design.
However, in conditions of limited availability of computational resources, and assuming that the model capacity will not compromise the solution accuracy (e.g. in the case of flow separation, pressure gradients, or rotation), more economical CFD analysis based on the sk − ε model with a calibrated value of mass transfer coefficient might be considered acceptable.

Conclusions
Work towards development of the unequivocal CFD modeling framework for AS tanks is in progress. The objective of the research presented in this paper was to select and experimentally validate a complete, transient numerical model that allows prediction of the hydrodynamics and mass transfer in a lab-scale aeration tank. The following key conclusions are drawn: • Considering two-phase aerated systems, accurate simulation of the liquid velocity field and turbulence parameters is of crucial importance for correct prediction of the air holdup and mass transfer coefficient. • The sk − ε and IAC models yielded the same air velocities and holdup as the model assuming fixed bubble size. The inaccurate prediction of inertia-dominated rotating flow and the use of an empirical equation for ε led to overestimation of the coalescence effects and mass transfer coefficient k L a 20 . • The SST k − ω + IAC approach accounting for coalescence and breakage effects was found to be the most accurate in reproducing the measured air velocity field and k L a 20 . • For the lowest operating airflow rates, the k L a 20 values obtained with the SST k − ω model and fixed bubble size were closer to the experimental results than the sk − ε model outcomes were. • The results obtained with the SST k − ω model highlighted the necessity of correct estimation of the bubble sizes, and the empirical adjustable model constants C and B , which, along with ε, play an important role in predicting the actual coalescence and breakup rates.
• Assessment of the oxygen transfer rates in AS process tanks is of crucial importance in the context of biochemical conversion reactions yield and optimization studies. In complex flow systems, such as those found in mechanically agitated aeration tanks, reliable determination of the k L a values relies on accurate prediction of the local hydrodynamics, accounting for the effects of turbulent interactions between the phases and resulting bubble sizes. Despite the high computational cost involved, it is recommended that the optimal modeling scheme for such scenarios is the SST k − ω model and IAC. In conditions of low operating air superficial velocities, the necessity for bubble coalescence and breakup models can be overcome through experimental characterization of the bubbly flow regime and resulting bubble size distribution in the tank. Where computational resources are limited, and when considering less complex flow systems (e.g. no swirl/rotation) a more economical CFD analysis based on the sk − ε model with calibrated mass transfer coefficient and bubble sizes may be appropriate. • The limitation of the work presented in this paper is its validity for lab-scale clean water tanks. Accordingly, follow-up work should aim to evaluate the applicability of the proposed IAC-SST k − ω modeling scheme for the AS tanks (full-or downscaled). Therefore, to correctly represent flow patterns induced by the aeration and mixing devices, the 'core' transient model may need to be expanded through implementation of the calibrated submodel, adjusting density in such a way as to account for the presence of MLSS.