Two-Phase Bubble Columns: A Comprehensive Review

We present a comprehensive literature review on the two-phase bubble column; in this review we deeply analyze the flow regimes, the flow regime transitions, the local and global fluid dynamics parameters, and the mass transfer phenomena. First, we discuss the flow regimes, the flow regime transitions, the local and global fluid dynamics parameters, and the mass transfer. We also discuss how the operating parameters (i.e., pressure, temperature, and gas and liquid flow rates), the operating modes (i.e., the co-current, the counter-current and the batch modes), the liquid and gas phase properties, and the design parameters (i.e., gas sparger design, column diameter and aspect ratio) influence the flow regime transitions and the fluid dynamics parameters. Secondly, we present the experimental techniques for studying the global and local fluid dynamic properties. Finally, we present the modeling approaches to study the global and local bubble column fluid dynamics, and we outline the major issues to be solved in future studies.


Introduction
Bubble columns are multiphase reactors where the gas phase is dispersed into a continuous phase (e.g., a liquid phase in two-phase bubble columns or a suspension in slurry bubble columns) in the form of a "non-coalescence-induced" bubble or of "coalescence-induced" structures. In particular, two-phase bubble columns are widely used in petrochemical, biochemical and chemical industries owing to their advantages in both operation and system design (viz. the simple system layout) [1,2]. The simplest bubble column layout consists of a vertical cylinder without internals, in which the gas enters at the bottom-through a gas sparger (i.e., porous gas sparger, perforated plate, ring gas sparger, needle gas sparger or spider gas sparger, ref. [3]) and the liquid phase is supplied in batch mode or it may be led in either counter-currently or co-currently to the upward gas stream. Eventually, reactions between gas and liquid components may be observed in some practical cases (viz., hydrogenation oxidation, phosgenation, chlorination, alkylation, . . . [4]). Despite the simple system design, bubble columns are characterized by complex fluid dynamic interactions [5]. For this reason, their correct design and operation rely on the precise knowledge of the fluid dynamic phenomena on different scales: (i) "molecular-scale", (ii) "bubble-scale", (iii) "reactors-scale", and (iv) "industrial-scale" (see Figure 1). At the "molecular-scale", fundamental chemistry is required to study the interfacial phenomena, catalysts and gas conversion processes, and to formulate mass transfer models (e.g., for example, the film theory, the penetration theory, the surface renewal theory, the boundary layer theory; the reader may refer, for example, to Figure 2 for a simple representation of these phenomena). At the "bubble-scale", experimental and theoretical works have proposed quantifying and understanding bubble size distributions (BSDs), bubble shapes, single bubble dynamics, "non-coalescence-induced" This review mainly focuses on the "bubble-scale" and the "reactor-scale" and aims to provide insights into the bubble column fluid dynamics, to support methods and criteria for scaling-up proposals. In addition, the "molecular-scale" is briefly discussed to provide insight and a brief overview of the mass transfer phenomena and modeling approaches. It is worth noting that the knowledge of fluid dynamics at the "bubble-scale" and at the "reactor-scale" can be quantified through the precise estimation of the local fluid dynamic properties (viz., the bubble aspect ratio, the bubble rising velocity, the bubble size distributions, . . . ) and the global fluid dynamic properties (viz., ε G , the mean bubble residence time, . . . ) fluid dynamic properties. In this respect, the correct quantification and estimation of these fluid dynamics parameters is very important for avoiding overestimations of investment costs (i.e., with respect to the downstream processing units and to the reactor itself) and design failures. The review is structured as follows: In Section 2, we discuss how the operating parameters, operating modes, liquid and gas phase properties, and design parameters influence the flow regime transitions and the fluid dynamics parameters. In Section 3, we present the experimental techniques to study the "bubble-scale" and the "reactor-scale" fluid dynamic properties. In Section 4, we present modeling approaches to study the "bubble-scale" and the "reactor-scale" bubble column fluid dynamics.

Bubble Column Fluid-Dynamics
In this section, the flow regimes, flow regime transitions, gas hold-up, bubble size distributions and shapes, local flow properties and mass transfer phenomena are discussed. We also discuss how the operating parameters (i.e., the operating pressure, p c , the operating temperature, T c , and the superficial velocities, U G and U L ), operating modes (i.e., the co-current, the counter-current and the batch modes), liquid and gas phase properties (i.e., liquid viscosity, µ L , the gas density, ρ G , the presence of active agents, . . . ), and design parameters (i.e., the gas sparger design, the column diameter, d c , and the aspect ratio, AR influence the flow regime transitions and the global and local fluid dynamics parameters(it is worth noting that the correct definition of AR relies on the initial liquid level (AR = H 0 /d c ), rather the column height, H c , as discussed by Sasaki et al. [6]). It is worth noting that the proposed discussion considers an air-water bubble column operating at ambient conditions as the "baseline" case; thus, the variations owing to the design and operation parameters are discussed starting from this the "baseline" case.

The Flow Regimes and the Flow Regime Tranistions in Bubble Columns
In the next sections, the prevailing flow regimes are described and subsequently, the relationships between the boundaries of the different flow regime transitions and the different parameters of the systems are described.

The Flow Regimes in Bubble Columns
The coupling between the gas and the liquid phases physically manifest in the flow regimes. Based on the previous literature as well as our experimental investigations, in the current authors' opinion, a change in the design parameters of bubble columns and a change in the phase properties affect the transition points between the flow regimes and do not influence the properties and main characteristics of the flow regimes themselves. Indeed, when the system design parameters or the phase properties are changed, the boundaries between the flow regimes change as a consequence. In particular, when considering the most general case, five flow regimes (and the flow regime transitions in-between) may be observed when increasing the gas flow rate (e.g., the superficial gas velocity, U G ) using a fixed system design and phase properties: (i) the mono-dispersed homogeneous flow regime; (ii) the pseudo-homogeneous flow regime; (iii) the transition flow regime; (iv) the heterogeneous flow regime; (v) the slug flow regime; and (vi) the annular flow regimes. The first flow regime observed is the "mono-dispersed homogeneous" flow regime, which is characterized by mono-dispersed BSDs, according to the change in the sign of the lift force coefficient (see ref. [7][8][9][10][11][12]); an example of this flow regime is the one observed by Mudde et al. [13], in which a "fine" gas sparger was used, or by Besagni et al. [14]. When increasing U G , the BSDs shift towards higher bubble diameters (according to the change in the sign of the lift force coefficient, and, thus, the "pseudo-homogeneous" homogeneous flow regime is observed; of course, if large bubbles are aerated [7] by using a "coarse" gas sparger, the lower boundary of the former flow regime moves towards a low U G and eventually, disappears. Generally speaking, the mono-dispersed and poly-dispersed BSDs are characterized by limited interactions between the bubbles and only "non-coalescence-induced" bubbles exist, as defined and observed by [7]. In general, one may refer to the homogeneous flow regime as the complete range of both mono-dispersed and pseudo-homogeneous flow regimes. With an increase in U G , the number of bubbles having a negative lift force coefficient increases and these continue migrating towards the center of the bubble column; at a certain point, the rate of this process leads to the formation of "coalescence-induced" bubbles [7] (as demonstrated by Besagni et al. [11,12]) when "coalescing" solutions are employed. Conversely, if "non-coalescing" solutions are employed, the transition flow regime is identified by the appearance of "clusters of bubbles", as defined in refs. [15,16] and observed by [7]. The flow regime transitions from the homogeneous flow regime towards the transition flow regime concept are displayed in Figure 3. Beinhauer [17] found that the first "coalescence induced" bubble appeared at about U G = 0.02 m/s and, by increasing the gas flow rate, both "non-coalescence-induced" and "coalescence induced" bubbles increased until the non-coalescence-induced" bubbles reached a maximum, and then decreased toward a constant value. Conversely, Scumbe and Grund [18], Krishna et al. [19,20] and Besagni et al. [7] reported that the non-coalescence-induced" bubbles increase until the flow regime transition, and then only "coalescence induced" bubbles increase. The difference is probably caused by the different gas sparger design ("fine-gas sparger" and "coarse-gas sparger"). At higher U G , the heterogeneous flow regime is achieved. This flow regime is associated with high coalescence and breakage rates and a wide variety of bubble sizes; in the heterogeneous flow regime, the liquid phase flows in an oscillatory manner, upwards and downwards, because the U G is not high enough to carry the liquid continuously upward. The reader should refer to the recent review proposed by Montoya et al. [21] for a complete discussion on this flow regime. When U G is increased, the slug flow regime may occur; however, in industrial applications, this flow regime is not observed, owing to Rayleigh-Taylor instabilities [22] which arise because of the large diameter effects. The quantification of Rayleigh-Taylor instabilities at the "reactor-scale" is quantified through the dimensionless diameter, D* H : In Equation (1), D H is the hydraulic diameter, d c is the bubble column (inner) diameter, σ is the surface tension, g is the acceleration due to gravity and ρ L − ρ G is the density difference between the observed that an increase in AR destabilized the homogeneous flow regime. Sarrafi et al. [29] observed no effect beyond 4 m. Of course, these results are related to the evalution of the BSDs inside the bubble columns [36]. Besagni et al. [37] found that, when using a "coarse-gas sparger", an increase in AR destabilized the homogeneous flow regime up to a critical value equal to approximately 5 in batch mode ( Figure 2) and 10 in counter-current mode ( Figure 3). Figure 4 summarizes the main results from the literature concerning the influence of the aspect ratio on the flow regime transitions; Figure 4 considers both (a) the regime transition between the pseuso-homoegenous flow regime and the transition flow regime and (b) the regime transition between the transition flow regime and the heterogeneous flow regime. Figure 4. Influence of the aspect ratio (AR) on flow regime transition points-data obtained by [28].

Gas sparger design.
It is known that the gas sparger design highly influences the flow regime transitions in bubble column. Generally, gas spargers with a hole diameter, d 0 , smaller than 1 mm (the "fine distributors", i.e., porous [32,33,38,39], membrane [40], needles [13,14,41] or sieve tray/perforated plates with small openings [37] see Figure 5 for a representation of these gas spargers) can maintain a stable mono-dispersed homogenous flow regime also at a high U G ; hence, the homogeneous flow regime is stabilized [13] because the "mono-dispersed homogeneous" flow regime exists [37] and the "pseudo-homogeneous" flow regime takes place after the former is destabilized. In contrast, using a "coarse gas sparger" (d 0 > 1 mm-see Figure 6 for a representation of these gas spargers) the "mono-dispersed homogeneous" flow regime may not exist and a "pseudo-homogeneous" flow regime may be observed at lower gas superficial velocities; finally, when using a "very coarse gas sparger" (d 0 >> 1 mm), the homogeneous flow regime may not be established and the transitional/heterogeneous flow regime may take place [42]. A well-designed "fine distributor" may stabilize the homogeneous flow regime far beyond expectations, as demonstrated by Mudde et al. [13]. An early study was conducted by Zahradnik [43], demonstrating that, in perforated plates, the smaller the gas sparger opening, the more the homogeneous flow regime is stabilized. Sarrafi et al. [29] found that the transition velocity decreases as d 0 increases, up to a value of d 0 = 0.0015 m. Thorat and Joshi [44] found that the homogeneous flow regime is stabilized with a decrease in d hole (d 0 = 0.0008-0.05 m). It is also worth noting that they observed that an increase in AR destabilized the homogeneous flow regime. Sal et al. [45] found that U G,trans decreases while increasing the gas sparger openings (d 0 = 0.001-0.003 m). Besagni et al. [46] compared two different "coarse-gas spargers" and found no significant differences between them; conversely, in more recent studies, Besagni et al. [14,37] compared the flow regime transitions in both "coarse-gas spargers" and "fine-gas spargers" (a needle gas sparger [14] or a perforated gas sparger [37]); they found that a uniform bubble bed stabilized the homogeneous flow regime, owing to the presence of the monodispersed homogenous flow regime.

Influence of the Bubble Column Operation Mode
Akita and Yoshida [47] (d c = 0.152 m, H c = 2.5 m) observed that the superficial liquid velocity (U L up to 0.04 m/s) had no influence in both the counter-current and the co-current modes. Jin et al. [48] (d c = 0.160 m, H c = 2.5 m) reported that the transition point is approximately the same for the three working modes if U L < 0.04 m/s, whereas, for higher U L (in co-current and counter-current modes), the transition velocity decreases when the U L increases. Otake et al. [49] (d c = 0.05 m, H c = 1.5 m) observed earlier flow regime transitions with an increased liquid flow rate in counter-current operations (U L up to −0.15 m/s). A similar conclusion was also drawn by Yamaguchi and Yamazaki [50] (d c = 0.04 m and 0.08 m). Extensive studies on the influence of the liquid velocity in the counter-current mode were performed by Besagni et al. [28,34,46,[51][52][53][54], who found that the counter-current mode destabilized the homogeneous regime in all of the different configurations studied (viz. the pipe sparger in the annular gap and in the open tube configuration [34] and the spider sparger [46]). Figures 7 and 8 summarize the main experimental results obtained by Besagni et al. [28,34,46,[51][52][53][54]. Recently, Trivedi et al. [55] performed a detailed study on the hydrodynamics of counter-current bubble columns and found results in agreement with Besagni and co-authors. It is worth noting that Besagni et al. [28] also studied the influence of the bubble column aspect ratio in the counter-current mode; the results of this study are presented in Figure 8, where it can be observed that an increase in AR destabilized the homogeneous flow regime up to a critical value, equal to approximately 5 in batch mode ( Figure 2) and 10 in counter-current mode ( Figure 3). Further explanation for this effect is provided in Section 2.2, when discussing the influence of the operation mode on the gas hold-up.  Influence of AR and U L on flow regime transition points-data obtained by [28].

Influence of the Liquid Phase Viscosity
The homogeneous flow regime may be either stabilized or destabilized depending on the liquid phase viscosity. This behavior was interpreted by Besagni and co-authors using the "dual effect of viscosity over the flow regime transition" [11] concept: (a) "moderate/high viscosities" stabilize the homogeneous flow regime, owing to increased coalescence [3,11,56] and the presence of cap-bubbles [11] (as experimentally observed by Wikinson et al. [56] and more recently, by Yang et al. [57]); (b) "low viscosities" stabilize the homogeneous flow regime, as the reduced coalescence increases the number of small bubbles; For example, Besagni et al. [11] (Figure 9a,b) found that µ L , depending on its value, either stabilizes or destabilizes the homogeneous flow regime compared to air-water systems (The Mono-Ethylene Glycol, MEG, concentration, c MEG = 0%: U G,trans = 0.0264 m/s, ε G,trans = 0.09; c MEG = 5%: U G,trans = 0.039 m/s, ε G,trans = 0.18; c MEG = 80%: U G,trans = 0.023 m/s, ε G,trans = 0.07). It is worth noting that the increased coalescence may also suppress the homogeneous flow regime and, for µ L > 8 mPa·s, it may not exist even with a 'fine gas sparger' [8,[58][59][60].
In this respect, Besagni et al. [14] observed that an increase in the liquid phase viscosity has an effect that can be approximated as the influence of an increase in the gas sparger opening-this conclusion suggests that these two effects contribute to a shift in the boundaries of the mono-dispersed homogenous flow regime to very low gas flow rates (Figure 9c,d). Ruzicka et al. [60] observed that the homogeneous flow regime is stabilized at low µ L values (µ L = 1-3 mPa·s) and that U G,trans increases with µ L in this range, whereas it decreases at moderate µ L values (µ L = 3-22 mPa·s). Olivieri et al. [61] observed a stabilization of the homogeneous flow regime up to µ L = 4.25 mPa·s, and then for higher µ L , a destabilization of the homogeneous flow regime. Finally, the above-mentioned results were also confirmed by ultrafast X-ray tomography investigations [62] (d c = 0.07 m, H c = 1.5 m, orifice gas sparger, d 0 = 1 mm). In this study, it was observed that, beyond µ L = 5.18-8.95 mPa·s, the homogeneous flow regime is no longer observed. It is worth noting that, in these systems, a non-Newtonian fluid behavior may exist; in this respect, we propose some speculations based on the conclusions of Olivieri et al. [61]; they interpreted the stabilization of the homogeneous flow regime as the relaxation time by considering the negative wake phenomena occurring in the rear of a bubble in the case of viscoelastic fluids. Please refer to Besagni et al. [11] for further discussion concerning this point. Figure 9. Influence of the liquid phase on flow regime transition points-ref. [12,14].
Influence of the Active Compounds: Inorganic Compounds ("positive surfactants") Inorganic compounds (i.e., electrolyte solutions) are attracted to the interface where they adsorb positively, lower the surface tension, cause coalescence suppression and, finally, stabilize the homogeneous flow regime [63][64][65][66][67][68]. Indeed, it is known that most electrolytes inhibit bubble coalescence in pure liquid phases [63][64][65][66][67][68] and, as a consequence, the homogeneous flow regime is stabilized [69,70]. For example, owing to their influence, the prevailing flow regime may change from heterogeneous to homogeneous, due to surfactant addition [71] as the boundary between the flow regimes changes. In this respect, the transition molar concentration, c t , is defined as the concentration, c, of the non-coalescent media above which bubble coalescence is drastically reduced. Depending on the concentration, n, we can define a "coalescent regime" (n* ≤ 1) and a "non-coalescent regime" (n* > 1), through the dimensionless concentration, n*-following the formulation of n t , proposed by Prince and Blanch [72] (a modified version of the formulation of Marrucci [63]): where B is the retarded Hamaker constant (B = 1.5 × 10 −28 J m), R g is the gas constant, T is the temperature of the system, r b is the bubble radius, and σ and ∂σ/∂n are the surface tension and the surface tension gradient with electrolyte concentration, respectively. This concentration has been found to be unique for each salt [73], valid for a swarm of bubbles and does not depend on U G [66,74]. Some of the relevant studies from the existing literature are here summarized, i.e., Besagni et al. [75] (d c = 0.24 m, H 0 = 3.0 m; kitchen quality NaCl- Figure 6a,b), Thorat and Joshi [44] (d c = 0.385 m, H c = 0.385-3.08 m, NaCl), Ribeiro and Mewes [76] (d c = 0.12 m, H c = 1.25 m, NaCl, Na 2 SO 4 , NaI), Kelkar et al. [77] (d c = 0.154 m, H C = 3.25 m), Grover et al. [78] (d c = 0.1 m, H C = 1.5 m. NaCl and CuCl 2 ), Rucizka et al. [70] (d c = 0.14 m, H 0 = 0.4 m; CaCl 2 ), Orvalho et al. [69] (d c = 0.14 m, H 0 = 0.4 m; Na 2 SO 4 , NaCl and kitchen quality NaCl). It is worth nothing that Rucizka et al. [70] and Orvalho et al. [69] observed a dual effect of the electrolyte concentration on the flow regime transitions. Figure 9a,b display the influence of NaCl concentration on the flow regime transition points; conversely, Figure 10a,b display the influence of AR and n* on the flow regime transition points, supporting the above-mentioned discussion on the influence of the concentration of the inorganic substances in the "coalescent regime" (n* ≤ 1) and the "non-coalescent regime" (n* > 1). Figure 10. Influence of AR and n* on flow regime transition points-Data obtained from ref. [28].
Influence of the Active Compounds: Organic Compounds ("negative surfactants") The organic substances (i.e., alcoholic solutions) are repelled from the bubble interface where they adsorb negatively, cause small increase of the surface tension and, finally, inhibit the coalescence [79,80]: the alcohol molecules are composed by hydrophobic and hydrophilic parts that are adsorbed at the interface when dissolved in water, causing the coalescence suppression [81]. Among the literature, some of the relevant studies are here summarized, i.e., Krishna et al [19,20,82] (d c = 0.1, 0.15 and 0.38 m, H c = 4 m, EtOH up to 1% vol ), Al-Oufi et al. [32] (d c = 0.102 m, H c = 2.25 m, EtOH up to 300 ppm wt ), Dargar et al. [83] (d c = 0.127 m, H c = 2.75 m, EtOH up to 5% wt ) and Besagni et al. [7] (d c = 0.24 m, H 0 = 3.0 m; 0.05% wt - Figure 6a,b). Figure 9a,b display the influence of EtOH concentration on the flow regime transition points; conversely, Figure 11a,b display the influence of AR and EtOH concentration on flow regime transition points. Figure 11. Influence of AR and n* (as defined in Equation (2) on flow regime transition points-data obtained from ref. [12].

Pressure and Temperature
An increase in pressure stabilizes the homogeneous flow regime [2,20,56,[84][85][86][87][88][89][90][91][92][93][94][95][96]. Indeed, increasing the pressure increases the break-up rate and reduces the coalescence rate: the new equilibrium between coalescence/break-up leads to a decrease in the bubble size and delays the appearance of large bubbles. The cause is the propagation of Kevin-Helmholtz instability and internal gas circulation, as stated in refs. [91,97]. In particular, Kitscha and Mustafaogullari [22] applied this instability to model the break-up phenomena in bubble columns operated at high pressures. It is worth noting that the Kelvin-Helmholtz instability can be classified as a surface instability and, thus, it is somewhat similar to the above-mentioned Rayleigh-Taylor instability, discussed in Equation (1). The main difference between these two instabilities concerns the relative velocities between the phases: the Kelvin-Helmholtz instability allows a relative velocity between the phases. Based on the Kelvin-Helmholtz instability, Wilkison and Van Dierendonck [98] proposed a maximum stable bubble size, in which all the disturbances in the liquid phase with a wavelength larger than the critical wavelength may break up a bubble. The interested reader may also refer to the discussion proposed by Lott et al. [99] concerning this instability. Also, the temperature seems to have a stabilizing effect over the homogeneous flow regime [78,100,101] due to the formation of small bubbles [1,2,101,102]. The combined effect of pressure and temperature was studied by Lin et al. [101] (p c = 0.1-15.2 MPa, T c = 298-350 K) showing a stabilization of the homogeneous flow regime; it can be stated that many industrial processes may be run in the pseudo-homogeneous flow regimes. Generally, at high pressures and temperatures and for p c > 10 MPa, the flow regime transition is mainly governed by the temperature, properties of the liquid phase, liquid velocity and gas sparger design [1].

The Flow Regime Maps
The flow maps help in predicting the prevailing flow regime based on the operating conditions and column geometry. An example is displayed in Figure 12. Lots of flow maps have been proposed for two-phase flows in large and small diameter pipes, and the most common map for a bubble column was proposed by Shah [24] (Figure 12a).  This map-limited to low viscosity systems in ambient conditions-presents the prevailing flow regimes (homogeneous, transitional flow and heterogeneous) based on the column diameter and superficial gas velocity. Conversely the flow maps produced by Zhang et al. [103], which also apply to air-water systems in ambient conditions, are limited to small-diameter columns but take into account both the gas and the liquid superficial velocities. The interested reader may also refer to the review by Wu et al. [104] for a comprehensive discussion on flow regime maps in vertical pipes and annuli. In the current authors' perspective, a general flow regime map is still missing; future studies should be devoted to proposing a comprehensive flow regime map based on the discussion of the flow regime transitions proposed in Section 2.2.1 and Figure 3.

The Gas Hold-Up
The gas hold-up, ε G , determines the residence time and, in combination with the BSD, the interfacial area (related to the reactor size); it is related to the column design, operation, phase properties and operating conditions. It is worth noting that, generally, the homogeneous flow regime is more sensitive than the heterogeneous flow regime to the operating parameters.

Influence of the Bubble Column Operation
Superficial Gas Velocity and Gas Holdup Curve Structure It is known that the gas hold-up curve (the relation between ε G and U G ) summarizes the complexity of the bubble column fluid dynamics, from the "bubble-scale" to the "reactor-scale" (see, for example, the discussion in ref. [69], Section 3.2). For example, "coarse gas spargers" lead to monotonic gas hold-up curves (Figure 13a), whereas, when using "fine gas spargers", a peak in the gas hold-up curve can appear (Figure 13b) owing to the mono-dispersed bubble size distribution, which leads to the hindrance effect (related to the Ledinegg instability) that is physically manifested by a peak on the gas hold-up curve, leading to a reversed S-shaped curve, as discussed and commented on in refs. [14,37] and displayed in Figure 3b.
Given the mass conservation to the gas phase, the gas hold-up is given by ε G = (U G /U swarm ); theoretically, if the bubbles travel at their terminal velocity, ε G will increase linearly with the gas flow rate. However, the coupling between the phases causes deviations from linearity (see ref. [60]). In the homogeneous flow regime, the hindrance reduces the bubble velocity, thus increasing ε G , whereas in the transition flow regime the "coalescence-induced" bubbles result in a decrease in ε G and cause ε G to decrease proportionally less than U G [34]. The presence of "non-coalescence induced" and "coalescence induced" bubbles has been investigated by different authors (i.e., [17][18][19]46]). For the sake of clarity, Figures Figure 15). Bubble columns may run batchwise (or, with U L < 0.01 m/s), in co-current or in counter-current mode [1][2][3]. Low liquid velocities generally do not affect ε G -as found by several investigators [24,47,[113][114][115][116][117][118][119]-because, if U L is low compared to the bubble rise velocities, the acceleration of the bubbles (caused by the non-stagnant operation) will be negligible [120]. For example, Akita and Yoshida [47] reported a negligible effect of liquid velocities up to 0.04 m/s, either in gas-liquid counter-current or con-current operations. Of course, at higher liquid velocities, the column operation influences ε G ; it is generally admitted that co-current operation tends to reduce ε G and counter-current operation increases ε G , as bubbles are either accelerated or decelerated by liquid motion [1,114]. Baawain et al. [121], showed that the counter-current or co-current operating mode influenced ε G for about 5% of its weight, and less than 1% of its bubble size, showing that the effect observed is mainly caused by the bubble rise velocity and not only caused by the bubble size. In agreement with this explanation concerning the influence of U L on bubble motion, different studies have reported ε G decreasing in co-current operations [49,90,[122][123][124][125][126][127][128] and increasing in counter-current operations [48,49,128]. Biń et al. [128] compared the three-operation mode showing that ε G increases with U L in counter-current mode and decreases or remains constant in co-current mode. The effect is more pronounced at a high gas velocity and the difference in ε G between co-current and counter-current modes is around 10%. The same trends were observed by Jin et al. [48], with a maximum difference of 2% between counter-counter and co-current operations and it was found that the influence of working mode is lower at high ε G . Similar trends were found by Otake et al. [49], as already discussed when considering the flow regime transition. Besagni et al. [34,46,[51][52][53][54] found that when increasing the liquid flowrate, a faster increase in the hold-up is observed at low U G , and the transition point also moves toward lower superficial gas velocities (Figures 16-18). This change was explained by the effect of the liquid flow, which slows down the rise of the bubbles, leading to higher hold-up-a more compact arrangement of bubbles leads to an earlier flow regime transition. Above a certain hold-up (depending on U L ), the liquid superficial velocity has no more influence on the hold-up; this is possibly caused by the limited influence of the operation mode on the fully established heterogeneous flow regime. In agreement with the above-mentioned study, Jin et al. [48] also observed that the influence of working mode is lower at high ε G .

Influence of the Bubble Column Design
Herein the influences of the main bubble column design criteria are discussed. The reader may also refer to the pioneering study of Wilkinson et al. [56] and to our recent paper concerning scale-up criteria [37].

Column Size
In a small diameter bubble column (Equation (1)), the wall effects alter the bubble size, rising velocity and liquid recirculation; therefore, when increasing d c , ε G may decrease [20,30,[129][130][131]. The increased recirculation was shown by Krishna et al. [130]. Urseanu et al. [131] observed the decrease in ε G (d c between 0.15 and 0.23 m) while working with viscous fluids, in agreement with Behkish et al. [132], who noticed that the effect of the column diameter on viscous fluid is higher. The data collected by Fai et al. [105] and Yoshida and Akita [112] show that the effect of the column diameter on ε G is negligible for columns larger than d c = 0.15 m. Hughmark [133] found an effect of column size on ε G up to a diameter of 0.10 m. Kata [134] conducted measurements in 0.066 m, 0.122 m and 0.214 m columns and found that ε G increases with a decreasing d c . Koide et al. [135] measured ε G in an 0.55 m column and found no significant difference from the literature values reported for columns less than d c = 0.60 m. Deckwer [136] found a difference in ε G between a d c = 0.041 m column and a d c = 0.10 m column. Hikita et al. [80] measured ε G in a d c = 0.10 m column and compared their results with the ones reported in the literature for columns with d c > 0.10 m, finding no appreciable effect of d c on ε G . Gopal et al. [106] measured the ε G in d c = 0.2 m, d c = 0.6 m and d c = 1.0 m columns and concluded that d c and the gas sparger do not significantly influence ε G . Nottemkamper et al. [137] measured ε G in d c = 0.19 m, d c = 0.45 m and d c = 1.0 m columns and obtained comparable results for the d c = 0.19 m and d c = 0.45 m columns but lower ε G values for the d c = 1 m column at high gas rates, which they attributed to its larger diameter. Koide et al. [138] observed smaller ε G values in columns smaller than d c = 0.2 m. Accordingly to the model proposed by Lemoine et al. [139], d c has an influence in columns of up to d c = 0.70 m in either the homogeneous flow regime or in the heterogeneous flow regime. Despite some contradictory results, a generally accepted rule of thumb is that d c = 0.15 m is large enough for the results to be scalable [1,24,98,140], as supported by different investigators [8,31,141]. This scale-up criteria relies on the fluid dynamics phenomena and the coupling between the gas and liquid phases, which can be interpreted by considering the instabilities previously described in Equation (1).

Aspect Ratio
The coalescence phenomena, the local fluid dynamics (see the local results described by Xue et al. [35]), the end-effects (i.e., top column effects and near gas sparger effects) all tend to destabilize the homogeneous flow regime and decrease ε G and are more evident in low AR bubble columns (see the experimental data summarized in . In particular, in systems where the bubbles are not at their maximum equilibrium size (and where coalescence may occur), the liquid height (H 0 ) will influence the extent of the coalescence. Hence, ε G decreases with H 0 , because the higher the column, the longer the time the bubbles have to coalesce and thus, the lower the mean residence time of the dispersed phase. In one of the very first studies, Yoshida and Akita [112] and Patil et al. [142] did not observe any significant effects of AR on ε G . Wilkinson et al. [56]-when presenting the scale-up criteria-discussed the results obtained by Kastanek et al. [143]: the influence of AR is negligible for H C > 1-3 m and with AR > 5. Zahradnik et al. [8] found that ε G decreases and the homogeneous regime is destabilized while increasing the initial liquid level up to a critical aspect ratio, AR Cr ; the authors concluded that their results support the assumption of a negligible influence of AR on ε G and flow regime transitions for AR > 5. These assumptions were also confirmed by Thorat et al. [111], who found a negligible influence of AR on ε G for AR > 5 (air-water) or AR > 8 (non-coalescing liquid phase). Sarrafi et al. [29] compared their experimental results with other experimental data and excluded any effect of the initial liquid level on the flow regime transition when H 0 > 3 m. Ruzicka et al. [31] found that an increase in liquid height decreases ε G and destabilizes the homogeneous regime up to critical values. Sasaki et al. [6] have found that an increase in liquid height destabilizes the homogeneous regime and decreases ε G (AR up to 5). The reader may also refer to the recent discussion by Besagni et al. [37], who compared values from a large experimental dataset to demonstrate that AR Cr = 5 ( Figure 19a) is ensured only for pure liquid phase systems with very large sparger openings operated in the batch mode (i.e., AR Cr increases to 10 in counter-current mode ( Figure 19b) and in pure liquid phase systems, with small gas sparger openings operated in batch mode ( Figure 21)). The relationship between AR and the gas hold-up in batch mode is displayed in Figure 19a; conversely, Figure 19b displays the relationship between AR and the gas hold-up in counter-current mode. Other data concerning the relationship between the gas hold-up and AR for different liquid phases are discussed in the next sections, concerning the influence of the active compounds. Figure 19. Influence of AR in large diameter bubble columns with spider sparger-data from ref. [28].

Gas Sparger
The gas sparger design influences not only the flow regime transition, but also the ε G value and, specifically, the ε G curve (the relationship between U G and ε G ). It is difficult to provide a general rule because of the many parameters involved (i.e., gas sparger design and operating conditions) and the controversial results presented in the literature. Despite some discrepancies, we can state that, when using a "fine distributor", the ε G curve increases linearly in the homogeneous flow regime, then reaches a peak and then rises again [18][19][20]32,33,82,129,145]; conversely, when using a "coarse distributor", the ε G curve grows continuously [19,49,111,112,131,146,147]. This was also proposed by Urseanu [19] and Deckwer [3]. The different behaviors are due to the bubble dynamics (i.e., bubble formation at the gas sparger and coalescence/break-up phenomena); in "coarse distributors" there is a continuous appearance of large bubbles [19], whereas, in "fine distributors" large bubbles start appearing after the flow transition [18]. In accordance with this discussion, "fine distributors" have higher ε G compared to "coarse distributors" because of their narrower BSD (mainly, a lower rise velocity) [8,129,148] and, for the same distributor, a decrease in d 0 leads to an increase in ε G (mainly, a lower number of large bubbles) [30,45,123,149]. In the current authors' opinion, these considerations can be explained by the flow regimes described in Section 2.2.1 and displayed in Figure 3.

Viscous Media
In the literature, both the increase and the decrease of ε G were observed as the liquid viscosity, µ L , increased. All the apparently contradictory results can be explained by interpreting them in terms of the "dual effect of viscosity over ε G , as described by Besagni et al. [11]-they found that ε G continuously (and non-linearly) increases when the MEG concentration increase, up to c MEG = 5% − µ L = 1.01 mPa·s, along with the contribution of small bubbles ( Figure 22a). Conversely, if the concentration is further increased from c MEG = 5% − µ L = 1.01 mPa·s to c MEG = 80% − µ L = 7.97 mPa·s, ε G decreases ( Figure 22b). For this last concentration, the ε G curve lies even below that obtained for pure water. The authors explained this behavior occurs because at low viscosities, the coalescence is limited, and the large drag force reduces the bubble rise velocity, causing an increase in ε G . When the viscosity increases, the tendency to coalescence prevails, creating large bubbles rising the column at a higher velocity, thus reducing ε G . In particular, the "dual effect of the viscosity" has been clearly observed in other experimental investigations [60][61][62][150][151][152][153]: Eissa and Schugerl [150] (d c = 0.12 m, H c = 3.9 m, d 0 = 2 mm) showed that ε G first increases (µ L < 3 mPa·s), then decreases (3 < µ L < 11 mPa·s), and finally becomes roughly constant (µ L > 11 mPa·s). Bach and Pilhofer [153] (d c = 0.10 m, d 0 = 0.5 mm) showed that ε G first increases (µ L < 1.5-2 mPa·s), then decreases (3 < µ L < 11-12 mPa·s), and finally becomes roughly constant (µ L > 12 mPa·s). Godbole et al. [152] (d c = 0.305 m, H c = 2.44 m, d 0 = 1.66 mm) showed that ε G first increases up to µ L = 2.23-4.75 mPa·s, then it decreases (7.81 < µ L < 52.29 mPa·s, depending on U G ), and finally becomes roughly constant (µ L > 52.29 mPa·s). Khare and Joshi [151] (d c = 0.20 m, H c = 3.0 m, d 0 = 2.0 mm) showed that ε G first increases up to µ L = 4 mPa·s, and then it decreases for 4 < µ L < 10 mPa·s. Ruzicka et al. [60] (d c = 0.14 m, H 0 = 0.2-0.8 m, d 0 = 0.5 mm) showed that ε G first increases for µ < 3 mPa·s, and it decreases for 3 < µ < 22 mPa·s. Olivieri et al. [61] (d c = 0.12 m, H c = 2 m, H 0 = 0.8 m, needle gas sparger, d 0 = 0.4 mm) showed that ε G first increases up to µ L = 4.25 mPa·s, and then it decreases at higher viscosities. Concerning the remaining literature surveyed, the main results are described here. Please note that these apparently contradictory results can be explained by the above-mentioned criterion. For example, Yoshida and Akita [47,112]  investigated the effect of µ L on the dynamics of bubble bed formation using several aqueous solutions of saccharose (1 < µ L < 110 mPa·s)-they found a decrease in ε G when the viscosity increased. Hwang and Cheng [154] (d c = 0.19 m, H c = 2.5 m, d 0 = 1 mm) investigated the ε G structure in highly viscous Newtonian and non-Newtonian media and observed that ε G decreases when the viscosity increases. Zahradnik et al. [8] (d c = 0.15 m, H 0 = 0.53 m, d 0 = 0.5 mm) found that moderate/high viscosities (3 < µ L < 110 mPa·s) decrease ε G . Yang et al. [57] (d c = 0.15 m, H c = 1.7 m, d 0 = 0.7 mm) investigated the influence of the viscosity (1 < µ L < 31.5 mPa·s) on ε G by using a viscosity-increasing agent; they observed that ε G decreases when the viscosity increases. It is worth noting that the influence of the viscosity may be also described in terms of non-Newtonian behavior (and vice-versa). This point has been discussed by Besagni et al. [11] by coupling the results obtained by Olivieri et al. [61] with the "dual effect of viscosity"; the non-Newtonian related stabilization of the homogeneous flow regime proposed by Olivieri et al. [61], along with the "dual effect of viscosity" relationship between the flow regime transition and the gas hold-up, help to explain the higher ε G for the non-Newtonian solutions observed by Godbole et al. [152].

Active Compounds
When considering active compounds, ε G increases mainly as a consequence of homogeneous flow regime stabilization for both electrolytes [8,76,77,79,80,[155][156][157][158] and ethanol [19,20,32,82,83,114,125,[159][160][161] solution. ε G also increases because the active material adsorbed at the surface is pushed towards the back of the bubble. This causes a surface tension gradient, opposing the tangential shear stress and thus increasing the drag and reducing the rise velocity. It is worth noting that some authors have also observed a dual effect for active compounds at high concentrations.

Inorganic Compounds
Previous literature focused on the influence of inorganic compounds has mainly concentrated on the relationship between ε G and the concentration ratio, c ≤ c t . Zahradník et al. [8] (d c = 0.14, 0.15 and 0.29 m, H c = 2.6 m) studied the influence of nine electrolytes and found that ε G grew continuously for c ≤ c t , but little change in ε G was observed for c > c t . These findings were also observed by Besagni et al. [28,162] in a large scale bubble column, operated both in batch mode and counter-current mode (see, for example, Figure 23). A similar dependence was found for all electrolytes at c = c t . Ribeiro and Mewes [76] (d c = 0.12 m, H c = 1.25 m) tested three electrolytes (NaCl, Na 2 SO 4 , NaI) with four concentrations in the "coalescent flow regime", and ε G was found to increase up to the critical concentration. Kelkar et al. [77] (d c = 0.154 m, H C = 3.25 m) reported an increase in ε G and a negligible effect of electrolytes (NaCl, CaCl 2 and Na 2 SO 4 ) above the critical concentration. Ruzicka and co-authors [69,70] (d c = 0.14 m, H 0 = 0.4, d 0 = 0.5 mm) observed a dual effect of inorganic compounds on ε G and flow regime transition with respect to the electrolyte concentration. The reader may refer to [75] as well as to refs. [69,70] for a further literature review. Figure 23. Influence of NaCl on gas hold-up (batch mode)-data from ref. [28].

Organic Compounds
Previous literature regarding the influence of ethanol has mostly focused on the relationship between ε G values and the transition point (see refs. [19,20,32,82,83])-regardless of the ε G value, all the studies agree in regard to the increase in ε G using ethanol solution. On the other hand, a limited number of studies considered the foaming phenomena which may occur when using organic compounds [161]. Indeed, foaming phenomena may occur as discussed in the experimental studies of Besagni et al. [7,12,37] and as predicted by the theoretical approach of Shah et al. [161]; in this respect, Figure 24 displays the gas hold-up curve in batch mode for the different aspect ratios (5 ≤ AR ≤ 12.5) and EtOH concentrations. Compared with the air-water system, the gas hold-up increases with the addition of ethanol; starting from a low ethanol concentration, the gas hold-up increases compared with the air-water case. If the ethanol concentration is further increased, foam phenomena are observed (i.e., Figures 24b and 25b-reference codes are in Table 3). It is also worth noting that foaming phenomena are not only related to the concentration, but to the aspect ratio itself. A discussion on these experimental results can be found in our previous studies. Figure 24. Influence of EtOH on gas hold-up-data from ref. [12]. Table 3. References codes for Figure 24 data.

R0
c EtOH,wt = 0%-Air-water gas hold-up curve R1 c EtOH,wt = 0.3%-Gas hold-up curve measured after waiting 30 s for every gas hold-up measurement point from low to high gas flow rate: run 1.

R2
c EtOH,wt = 0.3%-Gas hold-up curve measured after waiting 30 s for every gas hold-up measurement point from low to high gas flow rate: run 2.

R3
c EtOH,wt = 0.3%-Gas hold-up curve measured after waiting 120 s for every gas hold-up measurement point from high to low gas flow rate.

R4
c EtOH,wt = 0.3%-Gas holdup curve measured after waiting 120 s for every gas hold-up measurement point after each flow rate increase.

R5
c EtOH,wt = 0.3%-Gas hold-up curve when foaming phenomenon was observed: run 1 R6 c EtOH,wt = 0.3%-Gas hold-up curve when foaming phenomenon was observed: run 2 R7 c EtOH,wt = 0.3%-Gas hold-up curve when foaming phenomenon was observed: run 3 If the number of studies concerning foaming phenomena is limited, the number of studies concerning the influence of EtOh on gas hold-up are more common. Kelkar et al. [160] (d c = 0.154 and 0.3 m, H C = 2.44 and 3.35 m) studied the influence of ethanol (0.5-2.4% wt ) on ε G . Zahradnik et al. [159] (d c = 0.15 m, H C = 1.8 m) observed that the increase in ε G was higher with longer hydrophobic chains of n-alcohols. Rollbush et al. [114] (d c = 0.16 m, H C = 1.8 m) reported an increase in ε G using ethanol (1% vol ) up to 220%, compared with tap water. Hur et al. [163] investigated the effect of adding alcohol in bubble columns on ε G : The continuous mode decreased ε G in the homogeneous flow regime, whereas, in the heterogeneous flow regime, it increased ε G . The batch mode was less effective than the continuous mode in the heterogeneous flow regime. Pjoteng et al. [125] (d c = 0.10 m, H C = 1.8 m) reported an increase in ε G with the addition of ethanol (0.5% wt ), up to a pressure of 9 MPa. Guo et al. [164] studied a small diameter (d c = 0.1 m) and small-scale (H c = 1.8 m) bubble column. They found that ε G first increases, and then decreases, with an increasing ethanol concentration, owing to the non-linear effect of the alcoholic solution on the "coalescence-induced" and "non-coalescence-induced" bubbles. Similar results were presented by Syeda et al. [165] and by Shah et al. [161]. These studies [161,[164][165][166] explain the trend of ε G with respect to the mole fraction of one of the two components in the liquid mixture, based on the model proposed by Andrew [167].

Influence of Gas Properties
A variation in the gas properties may cause changes in bubble size distributions, bubble rise velocities and buoyancy [168]. Generally, the influence of the gas phase over ε G is discussed in terms of the gas density (ρ G ) [169]. A variation in the gas density may result from high pressure operation [2, 56,84,89,94,114,169,170], from the use of different gases [47,80,89,94,107,166,168,169,[171][172][173] or from a combination of both [89,94,169]. As the influence of the operating pressure is discussed in another section (Section 2.2.1), we here discuss the influence of the gas employed. Apart from Shulman and Molstadz [171]-who reported no effect of the gas properties-the literature agrees that, in columns of at least d c = 0.15 m, higher ρ G result in higher ε G [56]; a variation in ρ G means a variation in the residence time of the dispersed phase [114]. Also, the influence of ρ G over ε G has frequently been attributed to differences in bubble size distributions at the gas sparger [166,168,173]. Bagha et al. [166] (d c = 0.0382 m, H C = 1.14 m) and Koetsier et al. [172] reported that gases with higher density have higher ε G . Akita and Yoshida [47] (square column, 0.15 × 0.15 m) tested four gases (air, He, O 2 , CO 2 ) with water and showed that the effect of ρ G on ε G can be neglected, although ε G with helium are slightly lower for higher gas velocities. Hikita et al. [80] (d c = 0.10 m, H C = 1.50 m) observed an effect of gas density by testing different gases (air, H 2 , CO 2 , CH 4 , C 3 H 8 and H 2 -N 2 mixtures). Sada et al. [107] (d c = 0.073 m, H C = 0.95 m) concluded that the influence of ρ G on ε G should be considered in gas-molten salt systems because of the low gas-liquid density ratio at high temperatures. Özturk et al. [173] (d c = 0.085 m, H C = 0.95 m) considered different gases (Air, H 2 , He, N 2 , and CO 2 ) and observed that ε G increases when ρ G increases. Reilly et al. [94] (d c = 0.15 m, H C = 1.7 m) investigated various gases at atmospheric pressure (Air, Ar, He, N 2 , and CO 2 ) and conducted experiments at higher pressures (d c = 0.15 m) using N 2 and air in isoparaffinic solvents-higher ρ G extended the homogeneous flow regime and increased ε G , especially in the heterogeneous flow regime. Krishna et al. [89] (d c = 0.16 m, H C = 1.2-4.0 m) observed that ε G increases with ρ G as a result of both higher pressures and higher gas molecular weights (He, N 2 , Ar, CO 2 , and sulfur hexafluoride gases) in water. Jordan and Schumpe [169] (d c = 0.1 m, H C = 2.4 m) adjusted ρ G by increasing the pressure and by using mixtures of He and N 2 ; they concluded that the differences in ε G could be attributed only to differences in ρ G and not to the nature of the gas. Hech et al. [168] (d c = 0.15 m, H C = 2.2 m) tested different gases (air, He, N 2 , and CO 2 ) and found that that higher gas densities lead to higher ε G .

Pressure
The literature agrees that an increase in pressure results in higher ε G [1,20,56,84-86,88-90, 92-98,116,125,131,132,169,170,174-184]. Moreover, the influence of the pressure over gas hold-up is non-linear with U G ; whereas, in the heterogeneous flow regime, ε G always increases, in the homogeneous flow regime, some authors have observed an increase-even if less than in the heterogeneous flow regime-and others have observed no effect [1]. The literature also agrees when reporting a plateau for a certain value of pressure, above which, no significant effect of pressure on ε G is observed. The value of the plateau depends on the operating conditions [132], and typically, values are 6 MPa [185], 7 MPa [180] and 10 MPa [178]. The effect of pressure can be related to the liquid phase (i.e., µ L and ρ L ) and gas phase (i.e., ρ G ) parameters. The liquid phase has a limited influence [180]; the parameters of the gas phase have a higher influence [1]. It is worth noting that the gas density also has an influence over ε G at ambient pressure (Section 2.2.4). Ishiyama et al. [177] (CO 2 /water) reported a negative effect of pressure on ε G at pressures above 0.8 MPa in a heterogeneous flow regime; this behaviour was explained by an increase in µ L with dissolved CO 2 .

Temperature
Most of the authors reported a positive effect of the temperature, T c , over ε G [56,86,117,132,177,180,181,[186][187][188]; on the other hand, Pohorecki et al. [187] observed no effect of T c -possibly caused by the evaporation-and some authors observed a decrease of ε G T c increased [78,136,189,190]. It is worth noting that the studies reporting a decrease in ε G used at low U G and, in accordance with Yang et al. [190], turbulence is low under these conditions,. Under these conditions, an increase in T c , reduces the viscosity and weakly promotes turbulence; this promotes collision and increases film drainage speed and thus, coalescence. Grover et al. [78] (d c = 0.1 m, H c = 1.5 m, air-water) observed a negative effect of T c at atmospheric pressure (T c between 303 and 353 K), although for T c > 323 K, the effect of temperature is limited. Deckwer et al. [136] (d c = 0.041 m, N2/Paraffin system) observed a decrease in ε G until a plateau was reached. For the same range of U G and for the same d c as that of Deckwer et al. [136], Kölbel et al. [189] (d c = 0.041 m, H c = 4 m, H 2 , Ethylene/C 13-18 ) observed the same decrease in the homogeneous flow regime, but no effect in the heterogeneous flow regime. Yang et al. [190] (d c = 0.041 m, H 2 , CO/Paraffin) observed a decrease in ε G with an increase in both T c (293-523 K) and p c (1-3 MPa).

Gas Hold-Up Correlations
Different correlations have been proposed in the last few decades to correctly predict ε G in two-phase bubble columns. This section summarizes and described some frequently used ε G correlations. In particular, Table 4 summarizes the range of applicability of the reported correlations. In particular, in the following sub-sections, the main gas hold-up correlations are presented grouped in six main schemes of correlation: (a) the scheme of correlations by Lockett and Kirkpatrick; (b) the scheme of correlations based on ε G ; (c) the scheme of correlations by Akita and Yoshida [47]; (d) the scheme of correlations for Newtonian and non-Newtonian liquid phases; (d) the scheme of correlations based on the work of Syeda et al. [165]; (e) the scheme of correlations for bubble column scaling-up.

Scheme of Correlations by Lockett and Kirkpatrick
In 1975, Lockett and Kirkpatrick [191] proposed a fundamental scheme of correlations, which served as the basis for many subsequent studies; they tested different equations to describe ideal bubbly flow, and they found that for d b = 0.0005, none of the proposed equations are completely satisfactory. They concluded that the Richiardon and Zaki equation with exponent 2.39 is adequate if a correction factor is also included to take bubble deformation into account. Therefore, they proposed the following correlation: Scheme of Correlations Based on ε G : From Hugmark towards Reilly et al.
The first gas hold-up correlation presented was the well-known Hughmark [133] correlation; this correlation is based on his experimental data, which takes into account the effect of the liquid properties; this correlation was based on the observation that the gas hold-up for the air-water system in batch mode is correlated with the non-dimensional number, defined as U G [(1/ρ L )(72/σ)] 1/3 . The correlation reads as follows: Mashelkar [192] suggested a non-linear correlation to predict ε G values in air-water bubble columns in both the homogeneous flow regime and the transition/heterogeneous flow regime: As observed by [106], at low values of U G (in the homogeneous flow regime) the denominator is approximately constant and ε G increases linearly with U G . At relatively high values of U G (in the transition/heterogeneous flow regime), ε G decreases as U G increases. A scheme of correlation, similar to the one proposed by Mashelkar [192], was proposed by Kato and Nishiwaki [193], as reported by [24,193] and shown in Equation (6): where the terms in Equation (6) reads as follows: β = 4.5 − 3.5 − 2.548 · d c 1.8 and γ = 717 U G 1.8 /β .
Another modification of Equation (5) was proposed by Koide et al. [135] (Equation (7)), who conducted research to find the effect of the column diameter on flow properties using a commercial-scale, 5.50 m diameter column, approximately 7.00 m in liquid height. The measured values of the average gas hold-up show a certain amount of scatter, but when U G is less than approximately 5 [cm/s], they almost agree with the values calculated by Kato et al. [134] when d c = 5.50 m: where the terms in Equation (7) read as follows: Hikita and Kikukawa [194] correlated the measured ε G , in air-various liquid systems with Equation (8): Grover et al. (1986) [78] proposed a modified form of Hikita's correlation ( [80]) to incorporate the influence of temperature on the physical properties. Their correlation is given by Equation (9) where the constants, a and b, have been found to be 1.1 × 10 −4 and 5 × 10 −4 , respectively. The vapor pressure (P v ) of the solvent has been used as a correlating term to account for the temperature effect: Gestrich and Rähse [195] correlated the literature data for ε G and suggested the following equation, which takes into account the liquid properties, the column dimensions and the operating variables: (10) In Equation (10), the mean bubble diameter (d b ), usually ranging between 0.002-0.004 m, has been found to have no significant effect on ε G . Therefore, d b = 0.0003 m can be used in the equation.
Kumar et al. [196] presented the ε G data for air in several liquids and found that their own and the previous investigators' data can be correlated using Equation (11) as a function of the dimensionless gas velocity, as follows: As reported by [108], this expression should be used for gas rates below 0.10 [m/s]; at higher gas velocities it exhibits a maximum value.
Friedel et al. [197] reported a correlation for calculating the gas hold-up in downhole bubble columns. It has the form of Equation (12), as reported by [24]: where ε * = U G U G +U L . Hikita et al. [80] obtained experimental data on the fractional ε G in bubble columns with various gases and pure liquids or aqueous solutions of non-electrolytes as well as for air in aqueous solutions of various electrolytes. These data have been used to study the effects of the physical properties of gas and liquid on the gas hold-up. As a result, a new dimensionless correlation has been presented and reads as follows: where f represents the correction factor, which is the factor by which ε G is increased by the presence of electrolytes in water. It assumes different values depending on the ionic strength (I). As the exponents show, the physical properties of the gas phase are of little relevance, but, if omitted, the mean deviation of the correlated data rises to 15% as opposed to 4% when included. Such properties could affect the relative gas hold-up, especially under high pressure conditions, but no systematic research has yet been done on this particular aspect. Reilley et al. [108] used a statistical approach, instead of performing a dimensional analysis, to develop a generalized correlation for gas hold-up which would cover the entire range of the physical properties (u G , σ, ρ G , ρ L , µ G , µ L ) of both the gas and liquid involved in their study. The result of their analysis using a non-linear multi-regression technique is given by Equation (14), and is in good agreement with their experiment and also with thirteen literature datasets selected by the authors to correspond to its region of validity. For example, only data for columns of at least 0.15 m in diameter were considered: Elgozali et al. [198] fitted data on gas hold-up in the bubble column for all studied batches (coalescent and non-coalescent) and two types of ejector. In the non-coalescent system, the size of the primary gas bubbles generated in the ejector stayed almost the same during their journey along the bubble bed. In coalescent systems, the effect of primary gas dispersion was suppressed by fast coalescent at the column bottom. Therefore, they correlated the two subsets separately and they ended up with a power-law empirical formula (Equation (15)), where the kinematic surface tension The empirical constants are reported in Table 5. It can be inferred that the effect of liquid phase properties, characterized by the exponents b and c, respectively, is almost identical in both classes of batches, while a significant difference arises in exponent a for the superficial gas velocity, which depends on the flow regime. 4 : The Akita and Toshida Scheme Akita and Yoshida [47] measured the gas hold-up for different gas-liquid systems and by means of dimensional analysis, proposed the well-known gas hold-up correlation shown in Equation (6):

Scheme of Correlations Based on
In Equation (16), the c 1 = 0.20 for non-electrolyte systems, whereas c 1 = 0.25 has been suggested for electrolyte solutions which cause a larger ε G compared with non-electrolyte liquids (as discussed in the previous sections). Equation (6) is still considered the state of the art in gas hold-up correlations. Mersmann [199], as reported by [24], proposed a correlation similar to Equation (16) in [47]: A discussion concerning this scheme of correlation has been provided by Santus and Salvagno [200]. A similar scheme of correlation was used by Bach and Pilhofer [153], Riquarts and Pilhofer [201] and Iordache et al. [202]. In particular, Bach and Pilhofer [153], as reported by refs. [24,203], proposed a gas hold-up correlation for pure organic liquids over a wide viscosity range: Equation (18) does not apply to glycerol solutions because, as concluded by the authors themselves, in this system, the gas hold-up is substantially different since it exhibits a maximum at about 2-3 [mPa·s], owing to "the dual effect of the viscosity". A modified version of Equation (17) was proposed by Riquarts and Pilhofer [201], as reported by [204]; the proposed correlation for an air-water system reads as follows: Iordache et al. [202] interpreted results concerning the dependence of the bubble void fraction on the superficial velocity in a gas-liquid dispersion in a bubble column using a stochastic model for the flow phenomena, based onū • (viz. the mean individual velocity): Sada et al. (1984) [107] modified the equation used in [47] by including the term, gas-liquid density ratio, to better stress the influence of the gas phase: Scheme of Correlations for Newtonian and Non-Newtonian Liquid Phases Godbole et al. (1982) [152] determined the fractional gas hold-ups using the dynamic disengagement method in a 0.305 m diameter bubble column, and they presented empirical correlations based on data covering wide ranges of viscosities in Newtonian and pseudoplastic non-Newtonian solutions. From the measurements for water-glycerine solutions in the viscosity range 0.00423-0.246 Pa·s, they proposed the following correlation, which can be used to predict the the gas hold-up prediction in viscous Newtonian liquids: The gas hold-up data for CMC (Carboxymethyl cellulose) solutions have been correlated with another empirical correlation (Equation (22)), which is based on the data covering an apparent viscosity of 0.018-0.230 Pa·s. The apparent viscosity is given by µ L = K(5000.0U G ) n−1 , where K is the consistency index and n is the flow behavior index: Another correlation for highly viscous CMC solutions (≥0.02 Pa·s) has been proposed similarly to those presented by [205] for 0.102 and 0.14 m diameter bubble columns: Using the gas hold-up data for highly viscous CMC solutions (≥0.02 [Pa·s]) in the three bubble columns (the two columns used by [205] and the one used by [152]) Godbole et al. [152] proposed the following empirical correlation to also take into account the diameter effect: Viswanathan and Rao, on the basis of the cell model [206], obtained a correlation (Equation (26), reported by [204]), in which u br is the bubble rising velocity. It predicts a dependency of the gas hold-up on the column diameter which is similar to that in the model proposed by [204] but somewhat smaller: Kawase and Moo-Young [207] proposed the following empirical correlation (Equation (27)), where ν a is the apparent kinematic viscosity and the apparent dynamic viscosity is expressed in terms of the consistency index, K, as follows µ a = K . γ n−1 = K(2U G /d c ) n−1 : It is important to highlight that Kawase and Moo-Young [207] suggested that a more precise correlation should also account for the effect of the surface tension on the gas hold-up which they have not included. In addition, Kawase and Moo-Young [204] developed a theoretical model for gas hold-up in bubble columns with Newtonian and non-Newtonian fluids using the concept of a characteristic turbulent kinematic viscosity in bubble columns: Their correlation seems to be in reasonable agreement with the available data for Newtonian fluids over a wide range of column diameter (D = 0.1-1.07 [m]). However, to better correlate the experimental data for large bubble columns (D = 5.5 [m]) obtained by [208] which lie significantly above their model, Kawase and Moo-Young [204] proposed another correlation: Schumpe and Deckwer [203], based on a review of literature data as well as on their own experimental data, developed a single correlation for different viscous Newtonian and non-Newtonian solutions. Two hundred hold-ups were correlated using the same dimensionless numbers suggested by [47]. However, instead of the group ε G /(1 − ε G ) 4 the gas hold-up itself was empirically correlated for simplicity, leading to Equation (30): Equation (30) is suggested for the gas hold-up in the heterogeneous flow regime; conversely, in the slug flow regime, a theoretical approach has been found to apply, leading to Equation (31): Zou et al. [209] developed a gas hold-up correlation, implicating an effect of the operating temperature in the range of 25-96.56 • C. Following what was reported by [47,80], they neglected the effects of the nozzle diameter, the column diameter, the clear liquid height and the superficial liquid velocity. Thus, the factors affecting the gas hold-up were considered to be the superficial gas velocity, the liquid density, the liquid viscosity, the surface tension of the liquid, the gravitational constant and the liquid vapor pressure, which have been correlated by Equation (32): For the air-5% NaCl solution, the correction factor, f, has been defined as is the calculated values from Equation (32) for the electrolyte solutions. This correction factor has been found to depend on the operating temperature ( Table 6). The vapor pressure of liquid (P v ) was used to indicate the effect of temperature on the gas hold-up, since the change tendency of the gas hold-up with temperature is basically similar to the relationship between P v and temperature. An increase in temperature results in an increase in the gas hold-up, particularly above 75 • C, 65 • C and 80 • C for the air-water, air-alcohol and air-5% NaCl systems, respectively. Zou et al. [209] also observed that ε G (air-alcohol) > ε G (air-5%NaCl) > ε G (air-water). Finally, Kawase et al., in their review on bubble column reactors, [140] reported the correlation proposed by [210] and given by Equation (33): Scheme of Correlations of Syeda et al.
One of the most famous and promising schemes of correlation that also accounts for the "dual effects" was proposed by Syeda et al. [165]; they used the same four dimensionless groups proposed by [80] and adopted them into their correlation for pure liquids (Equation (34)) which includes the probability factor as the bubble stability criterion: The constant (c 1 ) and the exponent (b) have been determined experimentally in combination with the model proposed for binary mixtures (Equation (35)), where We 1 and We 2 represent the Weber numbers of the pure liquids with lower and higher surface tensions, respectively, in the binary mixtures: At x = 0 or 1, the correlation coincides with the model suggested for pure liquids (Equation (35)). By comparing the second equation with experimental gas hold-up data for four pure liquids and four organic binary mixtures, the values of c 1 and b have been found to be 1.334 and 0.032, respectively. The Weber numbers used in the two equations have been derived strictly for pure liquids. Changes in bubble size and bubble mobility that might affect the gas hold-up of the mixtures are included in the frothing parameter, crk 2 /σ.

Scheme of Correlations for Bubble Column Scaling-Up
In 1992, Wilkison et al. [56] proposed a well-known gas hold-up correlation for bubble column scaling-up. Later, Besagni and Inzoli [211] proposed a generalization of the original correlation which also takes the bubble column geometrical parameters into account. The scheme of the gas hold-up correlation reads as follows: where U "coalescence induced" bubbles and U "non-coalescence induced" bubbles and are the rising velocity of the "coalescence induced" bubbles and the "non-coalescence induced" bubbles, respectively. The different terms in Equation (36) were derived by changing the equation of Wilkinson et al. [56] as follows: In Equations (37)-(39), the non-dimensional groups are defined as follows: dimensionless diameter (D* H -Equation (1)), the dimensionless bubble column height (AR, aspect ratio refs. [28,144]), and the dimensionless sparger opening (d* o ); the last one is defined as the ratio between the bubble size produced at the gas sparger, estimated by the correlation in ref. [212], as the diameter of the gas sparger opening.  Hikita and Kikukawa, 1974 [194] (as reported by [24]).

The Bubble Size Distribution and Shapes
An understanding of the bubble properties, size distribution and shapes and rising velocity is of fundamental importance for understanding and modeling the flow dynamics and mass transfer phenomena in bubble columns (see, for example the review of Risso [213]). Indeed, in the current authors' opinion, the stabilization/destabilization of the homogeneous flow regime and thus, the changes in ε G are caused by the modifications to the BSDs connected to the bubble interface properties [65,214,215]. Besides this concept, using a practical point of view, apart from ε G , another fundamental parameter for the bubble column fluid dynamics is the bubble size distribution (BSD) which provides, along with ε G , an evaluation of the interfacial area [140] and is used in computational luid dynamics (CFD) for model set-up and validation [9]. The variation in the BSD is among the main reasons for the effects of the operating parameters on ε G and the flow regime transition (i.e., stabilization of the homogeneous flow regime, increase in ε G and, eventually, the plateau effect). In addition to the BSD, the bubble shapes should be considered (i.e., the aspect ratio). Indeed, the interface shape and size is important for characterizing multi-phase flows properly (i.e., heat and mass transfer at the interface and the closures in computational fluid dynamics). Generally speaking, the bubble shape depends on the system variables (as discussed by Haberman et al. [216]): the bubble rising velocity (u b ), the equivalent bubble diameter (d eq ), the difference between the phase densities (ρ L -ρ G ), the liquid viscosity (µ L ), the gas phase viscosity (µ G ), the surface tension (σ) and the gravity acceleration (g). These variables are coupled into non-dimensional numbers: • the Eötvös number: • the Morton number (defined only by the properties of the phases): • the Reynolds number: • the Weber number: These non-dimensional numbers are used for relating the bubble properties, as shown by Clift et al. [217] for bubbles rising in an infinite medium. Clift et al. [217] proposed a graphical correlation to relate the bubble properties in terms of Eötvös (Eo), Morton (Mo) and Reynold (Re b ) numbers. For a constant Morton number, the bubble shape evolves from spherical to ellipsoidal to cap-shaped when increasing its equivalent diameter (corresponding to the Eötvös number) (see, for example, the flow visualistaions in refs. [7,11]). The bubbles may be considered spheres when the surface tension and/or the viscous forces are larger than the inertial forces. Ellipsoidal bubbles are defined as oblate spheroid (an axisymmetric ellipse with respect to the minor axis). For example, Chao [218] argued that air bubbles in the water phase maintain a spherical shape up to Re = 400. When Re > 400, bubbles become flatted: first spheroidal and subsequently, spherical-cap bubbles. The rule of thumb is that bubbles having aspect ratio larger than 0.9 are spherical. Several attempts have been made in the literature to correlate the aspect ratio to dimensionless parameters (see refs. [219][220][221]): Eo [222,223], We [223][224][225], Tadaki number [217,226], the Mo and the Eo numbers [227] and the Re and the Eo numbers [228]. Unfortunately, most of these correlations were developed for single bubbles/drops, and they may not be suitable for dense bubbly flows, as discussed by Besagni et al. [11,53,229]. To this end, Besagni et al. [12,229] proposed a correlation for the bubble aspect ratio (E) for dense bubbly flows which reads as follows: In Equation (44), a is the major axis of the bubble; in particular, the coefficients in Equation (44) were obtained by experimental data and are listed in Table 7. Another correlation for dense bubbly flows, valid only for air-water systems, was proposed by Besagni and Inzoli [53]; it has been further validated by Besagni et al. [230] and reads as follows: The reader should refer to Besagni et al. [230] for a discussion concerning the prediction of the bubble aspect ratio in dense bubbly flows under different experimental setups. The precise understanding of the relationships between the bubble shape and size and the gas and liquid velocities provide valuable information concerning the relationships between the bubble shape and the "reactor-scale" flow conditions, as stated by Ziegenhein et al. [231]. Unfortunately, as discussed by Leoeonard et al. [1] the influence of the superficial gas velocity over d b is far from being understood and controversial results have been published. It is generally admitted that-in the homogeneous flow regime-d b increases [156] and, after the transition, "coalescence induced" bubbles appear, whereas the contribution of the "non-coalescence induced" bubbles remain constant [18]. Some authors reported an increase in d b with superficial gas velocity, both in the homogeneous and heterogeneous flow regimes [35,38,85,87,127,183,[232][233][234][235][236][237][238][239], while other authors reported no effect of U G over d b [188,240] and some others reported a decrease in the bubble size [241]. In the literature, both unimodal [239] and bimodal [242][243][244][245] BSD have been found depending on the gas sparger design and operating conditions. Lau et al. [242] and Wongsuchoto et al. [244] observed a transition from unimodal BSD to bimodal BSD with increasing superficial velocity of the air. This change in the distribution of bubble size has been justified by increased coalescence [242] or break-up [244]. Besagni and Inzoli [46,52,53,229] described the BSDs in a polydispersed homogeneous flow regime in the batch and counter-current modes. The changes in the BSD (from bimodal to unimodal) in the counter-current mode were based on the force pushing the small bubbles toward the center of the pipe. This has been confirmed by the numerical studies of Lu et al. [246] and Lu and Tryggvason [247,248]. Generally, lot of parameters influence the bubble size distribution, and it is difficult to provide a general rule for BSD prediction. Similarly, the influence of the superficial gas velocity over the bubble shape (viz. the bubble aspect ratio, E) is far from being understood. For example, Figures 25-27 present the experimental findings of the relationships between E and the gas and liquid superficial velocities, obtained by the PoliMi and HZDR research groups (see, for example, refs. [7,12,46,53,[229][230][231][249][250][251]). These experimental findings concern air-water flow; experimental results concerning other working fluids are summarized in ref. [12].   Generally speaking, from the above-reported experimental findings, it can be observed that a general relationship between Eo and E seems to exist; larger bubbles seem to be characterized by lower E and thus, are likely to be flattered (E is in the range of 0.4-0.7); conversely, the small bubbles have higher E and thus, are likely to be spherical (e.g., if d eq > 1 mm, E > 0.7). It can be observed that, in the different cases considered, changing U L and U G does not affect the E significantly. A detailed discussion of the above-mentioned experimental findings can be found in ref. [230] and thus, it is not repeated here. A valuable conclusion is that the relationship between bubble size and shape (in the boundaries of the homogeneous flow regime) mainly depends upon the liquid phase properties and not on the flow conditions (i.e., the gas and liquid flow rates). Generally, distinct differences are observed in smaller bubbles, whereas the aspect ratios approach each other with increasing bubble size. It is also worth noting that, in some flow conditions, it is found that as the bubble shape is not highly dependent on the flow conditions, it can be speculated that, in these cases, E mainly depends on the prevailing bubble size distribution, rather on the existing flow pattern [250].

Gas Sparger Design
The gas sparger design (i.e., the hole opening and the number of holes) along with the properties of the phases and the operating conditions determines the BSD at the gas sparger region that will, thus, evolve into the BSD in the developed region of the column. A review concerning the bubble formation in gas-liquid systems was proposed by Kulkarni and Joshi, to whom the interested reader should refer [252]. The role of the gas sparger in the BSD was investigated by Camarasa et al. [253] and Polli et al [254] for coarse and fine distributors and by Miyahara and Tanaka [255] and Koide et al. [256] for porous gas spargers. Generally speaking, a "coarse distributor" ensures a narrow/uniform sized BSD [8,13,129,148]. On the other hand, a larger opening would generate a larger bubble that would evolve into a log-normal or more complex BSD-see, for example, the brief discussion proposed by Besagni et al. [14]. An insight into the bubble formation process at the gas sparger is given by Hur et al. [257] for a single nozzle and by Kazakis et al. [258] for a porous gas sparger considering different liquid phases. Finally, Besagni et al. [7,46] investigated the BSDs near the gas sparger, the phenomena in the developing BSDs and compared their results with correlations from the literature.

Viscous Media
Generally speaking, high viscous media are characterized by more stable interfaces and thus, promote the formation of large bubbles at the gas sparger [3,259], bubble coalescence [3,56,58,59] and decreases in the bubble breakup rate [24,56]. Viscous liquid phases also change the nucleation process [260] as the bubble detachment in viscous liquid phase ins mainly determined by the viscous forces, rather than the inertial and surface tension forces [261]. In highly viscous media, a bimodal BSD is observed-for example, d eq < 1 mm and d eq > 20 mm [262], 1 mm < d eq < 10 mm and 10 < d eq < 150 mm [152], 0.7 mm < d eq < 10 mm and d eq > 10 mm [57]. Rahba et al. [62] reported bimodal BSD with large bubbles having d eq = 40-45 mm. For example, for viscosities beyond µ L > 30 mPa·s, the contribution of the small bubbles results in a further increase in the total ε G with increasing viscosity [59], instead of the levelling at a constant value, as reported by Eissa and Schugerl [150]. An exception to the above references is the study of La Rubia et al. [241], who observed a decrease of d eq from 4.6 to 4.2 mm while increasing U G . Besagni et al. reported a "dual effect of viscosity on BSDs"; the BSDs shifted towards low equivalent diameters for "low/moderate" viscosities and, for higher viscosities, cap-bubbles appeared in addition to the spherical-ellipsoidal bubbles. This behavior has been explained by two simultaneous effects: (i) at higher viscosities, coalescence is promoted and large bubbles appear; (ii) moderate viscosities stabilize and increase the boundary layer thickness between the bubbles, resulting in a decrease in the coalescence rate of the small bubbles (the momentum of the bubbles could not overcome the layer thickness). Similar behavior was observed by Yang et al. [57] in small scale bubble columns, using a "fine gas sparger"; they interpreted the increase in the number of small bubbles at moderate viscosities using the stabilization effect of the bubbly layers [60]. It is worth noting that increasing the liquid phase viscosity not only changes the prevailing bubble shape, but also the bubble dynamics (e.g., rising and horizontal velocities and bubble trajectories), owing to the balance between viscous and non-viscous forces [263]; the changes in the bubbly dynamics have been observed also from the point of view of aggregated bubbles [264]. Unfortunately, up to now, detailed experimental and numerical investigations on the bubble interface in viscous media have been quite limited, and there is a lack in fundamental knowledge, as reviewed in refs. [265,266]. Recently, Orvalho et al. [267] found that the bubble contact time increases with µ L , passing through three phases: (a) a rise at low µ L ; (b) a jump at intermediate µ L ; (c) a plateau at high µ L . This result seems to contradict with the above-mentioned findings and may be related to the diffuser flow conditions (i.e., the differences between bubble rising the bubble column and experiments on pair-wise bubbles); in regard to this concept, the interested reader may refer to the experimental study of Sun et al. [268]. Of course, future research studies should be devoted to performing 3D numerical simulations (e.g., see refs. [269]) and obtain local measurements, to study the bubble interface in viscous liquid phases and the bubble shapes and dynamics (e.g., see the pioneering study of Bhaga and Weber [166] and more recent papers [259,268]).

Active Compounds
When using active compounds, the bubble size distribution changes because of coalescence suppression, and a decrease in d b is expected (along with an increase in ε G ) [8,38,56,80,108,127,149,161,234,241,253,270,271]. Among the broader literature on the topic, we present some examples for the sake of clearness. Keitel et al. [270] found that the BSD shifts toward lower diameters when employing non-coalescing solutions. This was also observed by Shah et al. [161] (ethanol-water mixtures) and by Machon et al. [271], who observed a dependency on the active compound concentration. They found that, depending on the mixture concentration, the bubble diameter in water is much higher than in electrolyte and alcohol solutions (c > c t ) or, if c ≈ c t , the diameter is midway between the one in water and the one for c > c t .

Influence of the Gas Properties
It is thought that a possible cause for the differences in ε G , when changing the type of gas, is due to changes in bubble size distribution [166,168,173]. Koetsier et al. [172] postulated that argon coalesces less than helium near gas spargers. Özturk et al. [173] linked the higher ε G to smaller bubbles formed at the gas sparger. Hecht et al. [168] concluded that the most probable causes of the differences in ε G are the changes in bubble size distributions (caused by the changes in the break-up and coalescence rate) as well as dissimilar bubble rise velocities.

Influence of the Pressure and Temperature
Pressure An increase in pressure results in a decrease in bubble size (d b ) [1], and some authors have reported a decreased bubble velocity and an increased number of bubbles [87,98,272]. The multiple effects of pressure were also found by Jiang et al. [178], who observed a plateau for ε G at approximately 10 MPa and a plateau for bubble diameter at approximately 1.5 MPa; this suggests that pressure has another effect that decreases bubble diameter. Lin et al. [180] (up to 16 MPa) observed a decrease in  [91,92,170] (p c = 0.1-1.3 MPa) reported that the pressure has an influence over large bubbles only, suggesting that, in the homogeneous flow regime, the influence of pressure may be also negligible. The influence of pressure over d b is related to the gas sparger employed-Maalej et al. [182] (p c = 0.1-5.0 MPa) reported a lower influence of pressure for porous gas spargers than for perforated gas spargers (because of the narrow BSD in "fine distributors"). Kölbel et al. [189] (p c = 0.1-0.6 MPa) attributed the absence of a pressure effect in both the homogeneous or heterogeneous flow regimes to the narrow BSD at the gas sparger. Oyevaar et al. [185] (p c = 0.1-8.0 MPa) observed higher effects of pressure in perforated gas spargers than in porous gas spargers. Idogawa et al. [176] (p c = 0.1-15.0 MPa) found that the influence of the gas sparger becomes lower as the pressure increases. The effect of the pressure is also influenced by the other operating parameters, such as liquid phase viscosity, the chemical reactions and high temperature. In these cases, the effect of pressure may be reduced by the opposite influence of the other parameters. Urseanu et al. [131] (p c = 0.1-1.0 MPa) reported a minor effect of pressure on ε G while working with viscous fluids, because of the opposite effect of pressure and viscosity on d b . Ishibashi et al. [273] (p c = 16.8-18.6 MPa) observed no effect of pressure in the presence of chemical reactions. Considering previous high pressure and temperature studies, it seems that there is a negligible influence of pressure, as reported by Soong et al. [188] (T c = 293-538 K, p c = 0.1-1.36 MPa), Sangnimnuan et al. [115] (T c = 437-537 K, p c = 4.5-15.0 MPa) and Pohorecki et al. [187,240] (T c = 303-433 K, p c = 0.1-1.1 MPa). These results may be explained as occurring due to the saturation of gas through the evaporation of liquid, resulting in an increase in d b [56,187,240].

Temperature
Increasing the temperature reduced the surface tension, thus enhancing the bubble break-up and moving the BSD towards a lower diameter. Also, an increase in the temperature has also an influence in viscosity, which is reduced, thus promoting coalescence [1]. Therefore, the reduced surface tension and viscosity have opposite effects, but generally, a decrease in d eq is observed with an increase in temperature [180,183,188]. Lin et al. [180] observed a narrower BSD as temperature increased from 290 K to 351 K. Shafer et al. [183] (T c = 298-448 K) noticed a decrease in d b as long as p c was high compared to the saturation pressure and the evaporation was negligible. The role of evaporation was also discussed by Pohorecki et al [187], who found d b to not be dependent on T c (in the range T c = 303-433 K). Soong et al. [188] reported a decrease in d b when the temperature increased from 293 K (d b about 1.5 mm) to 538 K (d b about 0.4 mm). The reader should also refer to the section concerning viscous media for an insight into the influence of the liquid viscosity over d b .

Local Flow Properties
Local flow properties (i.e., the liquid and gas velocities, the local void fraction, etc.) provide an insight into the bubble column fluid dynamics, which is needed to characterize the mass transfer phenomena [274]. When considering the local void fraction, we refer to the local value of ε G , which is a function of both the radial and axial positions in the bubble column [275]. The axial and radial variation in ε G generate pressure changes, resulting in liquid flow recirculation. Therefore, the investigation of the local void fraction profiles provides an evaluation of the bubble column liquid phase recirculation [276]. The magnitude of the void fraction radial gradients and liquid velocity driven by the gas phase depend on U G , U L , the column design, the phase properties and the operating conditions. A review of the experimental data concerning the local void fraction profiles have been presented by Joshi et al. [276]. In general, three kind of profiles are found in the literature: wall peaked, flat and center peaked. In the homogeneous flow regime (and considering batch or co-current conditions), center peaked profiles occur because a large number of larger bubbles or agglomerates rise at a markedly higher velocity than smaller ones, and thus, they transport most of the gas, as reported by Besagni and Inzoli [34,46]. On the other hand, a wall peaked profile is due to a large number of small bubbles. This behavior is caused by the lift force, which changes its sign (for the air-water case, at d b = 5.8 mm), thus pushing the small bubbles towards the wall. This behavior is well known in the literature and was discussed by Tomiyama et al. [277] and Lucas et al. [278]. Increased gas flow rates resulted in greater profile curvature from the column wall towards the center as well as higher void fractions [34,46,125,276]. In particular, the existence of a pronounced radial ε G profile is among the main characteristics of the heterogeneous flow regime, which is characterized by strong liquid recirculation. Indeed, in this flow regime, the formation of "coalescence-induced" bubbles at higher gas flow rates leads to the above-mentioned increased curvature of the radial profiles. Conversely, flat void fraction profiles generally occur in the homogeneous flow regime using "fine gas spargers" under controlled conditions, as detailed by Mudde et al. [13]. In the literature, a large number of liquid circulation models have been proposed, for example, by Wu et al. [279], Nassos and Bankoff [280], Ueyama and Miyauchi [281] and Luo and Svendsen [282]. As stated before, the liquid velocity and liquid recirculation are strictly related to the local void fraction profiles. Also, the gas phase velocity is related to the local void fraction profile by means of the bubble size distribution.

Interfacial Area
The interfacial area (a i ) is one of the most important parameter for multi-phase reactors. Like the gas holdup and the BSD, it depends on the geometry, operating conditions and the phase properties. The specific interfacial area and ε G are related with: Besagni et al. [12], based on extensive experimental investigations of bubble sizes and gas hold-ups, obtained the results displayed in Figure 28. These data are a further confirmation of the above-discussed relationship between the liquid phase properties and the bubble column fluid dynamics. Figure 28. Influence of the liquid phase properties on the interfacial area-data from ref. [12].
In the literature, some correlations have been presented and have been compared by Besagni et al. [12]. The correlation proposed by Van Dierendonck et al. [283] reads as follow: where C = 2.5. In the case of electrolyte mixtures, C = 1.45. The correlation proposed by Gestrich and Krauss [284] takes into account the bubble column aspect ratio and reads as follows: The correlation proposed by Akita and Yoshida [285] takes into account the bubble column diameter and reads as follows: When ε G ≤ 0.14, Equation (49) is simplified as follows: Besagni et al. [12] proposed a modification of Equation (49) to take into account the bubble column aspect ratio. The proposed correlation, validated against a large set of experimental data, reads as follows:

The Mass Transfer
This section aims to discuss the underlying phenomena, theories and correlations to describe the mass transfer rate (k L a). It is worth noting that k L a is computed as the product of two parameters: (a) the liquid phase mass transfer coefficient k L ; and (b) the gas-liquid interfacial area (a i ). Accordingly, Martìn et al. [286] stated that the mass transfer rate is mainly determined by two mechanisms (viz. the bubble oscillations, related to the concentration profiles around the bubble, related to k L , and the contacting area, a). When considering previous literature, it is worth noting that bubble columns are mostly employed in slow reaction-absorption applications where the interphase mass transfer resistance on the gas side might be considered negligible compared with the gas-liquid side (k L ) [287]. Generally speaking, up to now, only a few theoretical correlations for k L a prediction in two-phase bubble columns [288], numerical models [289][290][291][292] or correlation-based predictions [139,293] have been used.

Physical Phenomena and Approaches
As previously mentioned (see Figure 2), the mass transfer phenomena occur across the bubble interface; thus, its correct estimatation relies on the precise understanding of the motion of bubbles in the liquid phase and the motion of the liquid phase around the bubbles. In this regard, the interfacial area and the shape of the interface are the fundamental parameters and details to correctly model. Generally speaking, interfacial turbulence promotes intensive mass transfer (in this respect, the interaction of micro-eddies and particles leads to an increase in turbulence intensity in the bulk phase [294]). The interested reader may refer, for example, to Lochiel and Calderbank [295] who proposed a discussion concerning the relationships between the mass transfer coefficient, the bubble shape, the steady state and the unsteady state phenomena. Different theories are available to explain the mass transfer rate (i.e., the film theory, the penetration theory; see the surveys in refs. [288,296,297]). In the film theory, the concentration distribution is linear and one-dimensional; the coefficient of mass transfer in this special case is proportional to the diffusion coefficient and inversely proportional to the thickness of the film. Conversely, the Higbie's penetration theory [298] assumes unsteady state absorption of the dispersed phase by a fluid element adjacent to the surface and moving at a uniform velocity during the mass transfer process. It is worth noting that, in this approach, an unhindered flow situation may be observed when considering the larger bubbles [299]. In this approach, the contact time is usually defined as the ratio of d b to bubble rise velocity; it characterizes the residence time of liquid elements at the interface and is used to estimate the liquid phase mass transfer coefficient (see, for example, refs. [300,301]); it is also worth mentioning that Nedeltchev et al. [288,296] assumed that the contact time depends on both the bubble surface and the rate of surface formation. Owing to the uncertainties in the estimation of the contact time, the penetration theory is inapplicable in different practical applications (e.g., turbulence may determine the renewal of the fluid layers surrounding the bubbles). It is also worth noting that Kawase et al. [302] proposed a semi-theoretical approach based on Higbie's penetration and Kolmogoroff's turbulence theories. Recently, Nedeltchev [297] proposed a modified version of the penetration theory that can be used to correctly predict the mass transfer coefficient in both the homogeneous flow regime and the transition flow regime.

Influence of the Bubble Column Design
Gas sparger opening. It is not surprising that the gas sparger has a large influence on k L a [24]; indeed, depending on the gas sparger openings, the prevailing flow regime changes (as previously discussed and observed by Besagni et al. [14,37]). The influence of k L a on the gas sparger openings has been observed, for example, by Koide et al. [138] and is displayed in Figure 29. k L a increases with a decrease in d 0 and increases (less than linearly) with an increase in ε G . The increase in the mass transfer with a decrease in the gas sparger opening is related to the prevailing flow regime determined by the gas sparger [14,37,45]-decreasing d 0 moves the prevailing BSD towards a lower bubble equivalent diameter and thus, increases the interfacial area. Hence, a good distribution of bubbles across the reactor increases the efficiency of the gas flow rate on k L a [286,303]. The relationship between k L a and ε G is related to the shape of the gas hold-up curve, depending on the flow regime; indeed, the bubble column fluid dynamics are highly influenced by the sparger design in the homogeneous flow regime, but not in the heterogeneous flow regime [28,111]. From these observations, it is clear that highly efficient bubble columns should be operated in the homogeneous flow regime, rather than in the heterogeneous flow regime, where the dispersed phase is transported through the column as large "coalescence-induced" bubbles. These "coalescence-induced" bubbles have a lower residence time, leading to a decreased conversion process. In other words, the further enhancement of ε G in the heterogeneous flow regime is marginal and certainly not commensurate with the imposed increase in gas flow rate. Figure 29. Influence of the gas sparger openings on k L a-data taken from Koide et al. [138].
Column diameter. Deckwer and Schumpe [287] found that the dependence of k L a on the column diameter exists only for column diameters smaller than 0.6 m. Conversely, Akita and Yoshida [47] ( Figures 30 and 31) found that the mass transfer coefficient is related to the bubble column diameter up to a diameter of approximately 0.15 m, which is in agreement with the "large diameter" concept described in Equation (1) and the experimental results of Vandu and Krishna [304] (Figure 32).

Influence of the Bubble Column Operating Conditions
In this section, the relationships between the operating parameters and k L a are discussed. In particular, the gas flow rate and the operating pressure are considered and commented on. It is known that the relationship between k L a and U G is well described by Equation (52): where n can be between 0.7 and 0.92 [285,302,305,306]. Thus, the relationship between k L a and U G is less than linear, in agreement with the considerations proposed in Section 2.6.1-the relationship between k L a and ε G is related to the shape of the gas hold-up curve, depending on the flow regime; indeed, in the transition/heterogeneous flow regimes, the dispersed phase is transported through the column as large "coalescence-induced" bubbles, having a lower residence time and thus, a lower mass transfer rate. This observation is in agreement with the experimental observations reported by Letzel et al. [170], Kang et al. [179] and Ozturk et al. [173] (Figures 33-35). On the other hand, it has been observed that increasing the operating pressure increases k L a, owing to the decrease in bubble size (as previously stated). Figure 36 proposes a summary of previous literature concerning the relationship between k L a and U G .     [47,138,173,304,307].

Influence of the Liquid and Gas Phase Properties
As previously stated, the liquid phase viscosity plays an important role in k L a. Indeed, the liquid viscosity, if above the maximum value dictated by the "dual effect", increases the mean diameter of the bubbles in the dispersion, thus reducing k L a. (as observed by Kang et al. [179]). In this condition, the lather bubbles are stable in the flow at higher viscosities. In addition, the liquid viscosity also decreases the liquid diffusivity (D L ) [299]. In the case of oxygen transfer, D L is proportional to µ L −0.57 [173].
It is worth noting that some authors have also considered the influence of the gas viscosity proposed by [306]. Finally, Ozturk et al. [173] (Figures 37 and 38) and Gopal and Sharma [106] (Figure 39) studied the influence of the gas phase on k L a and observed a lower influence compared with the above-mentioned operating parameters.

The Mass Transfer Correlations
In the last few decades, many correlations for mass transfer prediction have been proposed (see the literature survey proposed in ref. [308] and the discussion in refs. [309][310][311]); for example, the correlation of Schumpe [18] reads as follows: where K = 0.063 (water/salt solution), K = 0.042 (water, 0.8M Na 2 SO 4 ).
Most of the correlations [24,47,285] for k L a prediction in gas-liquid bubble columns neglect the effect of the bubble column aspect ratio; in this regard, Zhao et al. [312] proposed empirical correlations that consider the effect of liquid height on k L a. The correlations proposed by Akita and Yoshida [47,285] (Equation (54)) and by Shah et al. [24] (Equation (55)) read as follows: where D G,L is the diffusivity coefficient between the gas/liquid phases. In their empirical correlation, Öztürk et al. [173] used the surface-to-volume mean bubble diameter as a characteristic length rather than the column diameter: Finally, for the sake of clarity, other correlations are listed below: • Kawase and Moo-Young, 1987 [302]: • Hikita et al., 1981 [306]: • Kang et al., 1999 [179]:

The Experimental Techniques
A comprehensive review on the experimental techniques for multi-phase reactors was provided by Boyer et al. [313] and by Jakobsen [314], to whom the interested reader may refer for all the details concerning the broader framework of the measurement techniques for multi-phase flows. In this section, we provide an overview of the measurement techniques for both the local-scale and global-scale fluid dynamic parameters discussed in the previous section.

The Flow Regimes
It is well known that the correct definition of an approximate transition point helps in understanding and modeling the fluid dynamics behavior in multi-phase reactors [89]. Different methods have been presented over the last few decades, as discussed by Shaikh and Al-Dahhan [102], Nedeltchev and Shaikh [26] and Nedeltchev [27]. The most common methods employed are statistical methods (by interpretation of the ε G -U G curve), the swarm velocity [7,11,34,46,[51][52][53]75,76,89,91,229,315,316] and the drift flux [7,11,34,38,46,51,53,75,76,229,[316][317][318] approaches. These methods are based on the analysis of ε G data. Other flow regime identification methods are the fractal [319], wavelet [320,321], Hurst's [322,323], spectral and nonlinear chaos [93,320,[324][325][326][327][328][329][330][331][332][333] and Kolmogorov entropy [332][333][334] approaches. Some authors have used methods based on the Shannon entropy, which measures the degree of indeterminacy in a system [335,336]. Nedeltchev identified the flow regimes based on the information entropy theory [26]. Finally, the detection of the flow regimes using local measurements was demonstrated by Shiea et al. [337]. The interested reader should refer to the recent study proposed by Tchowa Medjiade et al. [334] (viz. power spectral density, standard deviation, fractal analysis, Kolmogorov entropy); they found systematic differences between the various methods. In the current authors' opinion, future studies should be devoted to proposing a comprehensive approach that is able to take into account all the possible flow regimes, without the possibility of subjective evaluation.

The Gas Hold-Up
The most common method, for measuring ε G is the bed expansion technique [3]. This method is based on the measurement of the liquid free surface before and after aeration. Another well-known method is based on the measurement of the differential pressure between two or more points in the column [80,86,114,116,122,123,125,126,135,141,170,179,194,198,209], under the hypothesis of negligible acceleration and friction pressure losses, as discussed by Leonard et al. [1]. This method has been applied by a large number of authors. For example, Hikita et al. [80] determined ε G by measuring the static pressure at three points in the column at 0.25 m intervals, the lowest being 0.15 m above the gas sparger. Another method is the dynamic gas disengagement proposed by Schumpe et al. [18], which may also be used for investigating the structure of ε G (in terms of "coalescence induced" and "non-coalescence induced" bubbles) and for investigating the flow regime transition.

The Bubble Size Distribution and Shape
Different non-intrusive and intrusive techniques have been proposed for measuring bubble shape and bubble size distribution [338][339][340][341][342][343][344][345][346][347][348][349]. Non-intrusive measurement techniques do not influence the fluid dynamics and for this reason, image analysis has attracted growing attention in recent years [350]. In the literature, image analysis has been applied to small-scale [242,351] and medium-scale bubble columns [183]; conversely, there is a lack of studies concerning bubble size distribution in large diameter bubble columns and large scale multi-phase reactors. Indeed, image analysis techniques are still limited in the resolution of large void fractions, highly unsteady flows and large bubble clusters [352]. Additionally, there are well-known issues associated with overlapping-if ε G exceeds 1%, more than 40% of the bubbles are overlapping [353,354]. Despite the proposed methods to deal with overlapping bubbles by different studies [242], no agreement has been found so far (i.e., some approaches may cause a reduction in the bubble sample size and/or an underestimation of the bubble diameter). Furthermore, by using the classical image analysis methods, only projected bubble images can be obtained; despite some proposals against this issue [355,356], a solution is far from being reached. For this reason, the use of 2D projected images is a common way of analyzing the bubble images [242,244,352,[357][358][359]. Concerning the approximation of the ellipse of the bubbles, Lage and Espósito [358] have stated that the error in the measurement of each axis of the ellipse is approximately 6%. In addition, if considering the error introduced by the hypothesis of oblate spheroid and the optical distortion, Lage and Espósito have estimated that the experimental error in the determination of d eq is in the range of 10-15%. Finally, when investigating bubble size distributions, the number of bubbles that need to be sampled to achieve a reliable BSD should be considered [360]. Various studies have sampled different numbers of bubbles-in the range of 50 and 300 [38,244,358,[361][362][363], and a sensitivity analysis on this issue has been presented in our previous paper [53].

The Local Flow Properties
Wire mesh sensors and the tomography of gamma rays, X-rays, electric or ultrasonic waves have also been used in previous literature. Boyer et al. [313] and Jakobsen [314] proposed a complete description and detailed discussion of these techniques. Hernandez-Alvarado et al. [364] compared many experimental techniques to evaluate the local flow properties in bubble columns; they found that pressure transducers and gamma densitometry were able to produce reliable measurements, regardless of the prevailing flow conditions. Another comparison between experimental techniques was presented by Singh et al. [365]. Optical probes were found to underestimate local void fractions. Conversely, electrical resistance tomography was found to be the least reliable method. Finally, the experimental findings regarding the wire mesh method were found to be highly related to the bubble size and mesh spacing. Generally speaking, needle, heat transfer and ultrasound probes are among the most common methods applied for studying the local flow properties in bubble columns [313]. Heat transfer probes are capable of simultaneously measuring the local void fraction and liquid velocity. Ultrasound probes are capable of measuring the local void fraction, bubble diameter and velocity. Needle probes are capable of simultaneously measuring the local void fraction, bubble chord lengths and rise velocities. In particular, two types of needle probes have been previously used for measurements in bubble columns: optical fibers and impedance/conductive probes. Optical fibers and impedance probes operate based on the differences in the refractive index or conductivity, respectively, of the liquid and gas phases. Signal fluctuations due to phase changes at the probe tip allow the measurement of local properties; these devices were developed to measure the rise velocities along with the other local properties [337,[365][366][367][368][369][370][371]. When considering the experimental error of the needle probes, two parts should be considered: (a) the statistical error associated with the measuring time (proportional to the square root of the measuring time); (b) the bias error due to the difficulties of piercing a bubble at the bubble edge [372] (viz. the blinding effect, drifting effect and the crawling effect). A complete discussion of these issues has been proposed in our previous paper, to which the reader may refer to, as well as previous literature (i.e., see refs. [127,355,[373][374][375][376][377][378]). Besides the local void fraction and rise velocities, probes provide chord length distributions (CLDs). Several authors have proposed methods to relate CLD to corresponding BSD. A process that deduces BSD from CLD is called "backward transform". In contrast, the deduction of a CLD generated from an existing BSD is called "forward transform". Liu and Clark [379] (assuming ellipsoidal and vertically rising bubbles) proposed a "backward transform" using an analytical relationship between probability density functions for CLD and BSD. Clark and Turton [380] proposed a "forward transform" to infer the CLD from the BSD for different bubble shapes. Hu et al. [381] used a "forward transform" to obtain the CLD from the log-normal and uniform BSDs. Al-Oufi [363] proposed an algorithm based on a double transformation (a backward and a forward transformation). This approach has been also detailed and applied in the dissertation of Carrasa [382]. Rudisuli [383] proposed a Monte Carlo method that simulates the bubble sampling process in the cross-section of the pipe. Hoang et al. [384] decomposed the CLD into sub-distributions which are considered to be the CLDs of bubble groups. These CLDs are converted into BSDs using the approach of Liu and Clark [379], and finally, the whole system BSD is constructed. It is important to notice that to relate the chord length distributions to the bubble size distributions, a correlation for the bubble shape is required. Therefore, a correlation between the aspect ratio (φ) and a bubble dimension parameter should be provided. Several attempts have been made in the literature to correlate the aspect ratio to dimensionless parameters (i.e., the Eötvös number, related to the bubble size). However, these correlations are obtained for single rising bubbles (or single drops in liquids [223]) and their use for bubbly flows is questionable. This issue has been deeply considered by Besagni et al. [229], to whom the reader should refer: the use of correct aspect ratio correlations is essential.

The Mass Transfer
The mass transfer coefficient, as reported by refs. [107,285]), can be determined from experiments involving absorption of a gas (e.g., oxygen or carbon dioxide) into various liquids. Using an experimental column that is operated continuously with respect to the gas flow, and batch-to-batch with respect to liquids, it is possible to write the following mass balance, on the assumption of complete liquid mixing, where *, as a subscript, stands for the equilibrium concentration:

The Modeling Approaches
Numerical modeling of bubble column reactors is a way of predicting the multiphase flow for improving the reactor design and understanding the reactor fluid dynamics. Many phenomenological models have been presented in the literature [281,309,[385][386][387][388][389][390][391][392][393][394], but they are limited in providing detailed information regarding the two-phase flow phenomena. On the other hand, computational fluid dynamics (CFD) is a promising method for studying the local and global bubble column fluid dynamics. Among the available modeling techniques, the Eulerian multi-fluid approach is the most common approach for simulating bubble column reactors, as reported in the review of Jakobsen et al. [395]. Different studies have applied the Eulerian-Lagrangian approach [396][397][398][399][400][401] but have still been limited to the simulation of small-scale reactors with low ε G . It is worth noting that, at present, it is possible to model the flow regimes listed in Section 2.2.1 by using different modeling approaches (see the discussion below); unfortunately, a comprehensive model, able to simulate the whole range of operating conditions in bubble columns, is still missing.

The Eulerian Multi-Fluid Approach
The Eulerian multi-fluid approach treats each different phase as interpenetrating continua and relies on an ensemble averaging of the multiphase Navier-Stokes equations [402]. Considering an isothermal flow without mass transfer, the URANS governing equations for the k-th phase are as follows: The terms on the right-hand side of Equation (9) represent the pressure gradient, the viscous and Reynold stresses, the body forces and the interfacial momentum exchanges between the phases. The above-mentioned model considers non-reactive systems; recently, other authors have proposed modeling approaches for reactive multi-phase flows in bubble columns, to whom the reader should refer [290][291][292]. The last term in the momentum equation comprises several independent physical mechanisms, as follows: The terms on the right-hand side of Equation (12) represent the drag, lift, virtual mass, turbulent dispersion and wall lubrication forces. The reader may also refer to the work of Vik et al. [292] for an interesting discussion concerning these interphase forces.
The accuracy of the simulations based on this approach is highly dependent on the closure models and on the numerical aspects (i.e., numerical methods, the time step size and the mesh). Following an overview of the advancements and current trends in bubble column simulation using the Eulerian multi-fluid approach is provided. The pioneering work, performed among in 1989 and 1999, reproduced the fluid dynamics in rectangular bubble columns using two-dimensional Eulerian simulations [403][404][405][406][407][408][409][410][411][412]. However, other authors [413][414][415][416] remarked that, for more accurate predictions of the bubble plume oscillation period and the turbulent quantities, three-dimensional simulations are required for both rectangular [413][414][415] and circular [416] columns. In addition, high order discretization schemes should be applied, regardless of the closure models implemented [417][418][419][420][421]. For example, Oey et al. [421] reported that numerical diffusion may prevent the transient nature of the two-phase flow emerging. In addition, beside the accuracy of discretization schemes, mesh resolution is a fundamental parameter in determining the capability prediction of the approach-despite the number of mesh sensitivity studies performed in the past, the proper mesh size to adopt in a bubble column reactor simulation is still an open debate [422][423][424]. It is worth noting that the mesh size also depends on the turbulence modeling approach implemented. Large eddy simulations (LES, [425]) of bubbly flows require a mesh size within both a higher and a lower bound in order to get a proper filter cut-off [426,427]. Conversely, RANS (Reynolds-averaged Navier-Stokes) simulations are less restrictive on the mesh size [10,414,420,[428][429][430]. Besides the numerical aspects, proper turbulence modeling is important for the accuracy of the simulations; multi-phase turbulence models are derived from their single-phase equivalents and terms modeling inter-phase interactions are added to the transport equations of the turbulence model (see refs. [413,431]). In this respect, the multi-phase equivalent of the standard k-ε model is widely used [413][414][415]428,[432][433][434][435]. Comparisons between the k-ε model and LES simulations have shown that similar performances are obtained in terms of average quantities; conversely, more accurate fluctuating quantities are achievable through LES [436][437][438][439], owing to the isotropic turbulence hypothesis of the k-ε model. The SST k-ω model has proven to be slightly superior to the k-ε model for simulating upward bubbly flow [440,441] and vertical pipes [10,430,[442][443][444][445][446][447]. It is worth noting that suitable prediction and selection of inlet turbulent boundary conditions are needed for reliable predictions [448]. In this framework, models based on the RSM (Reynolds Stress Models) approach seem promising [449][450][451][452] for improving the numerical prediction of the turbulence quantities to be used as inputs for bubble coalescencence and break-up models. In addition, suitable models that take into account the bubble induced contribution should be selected and developed [453][454][455]; for example, recent studies aimed to calibrate this contribution based on DNS (Direct numerical simulations) studies [456]. As previously mentioned, correlations for interfacial forces are implemented to model the inter-phase momentum exchanges (i.e., see the screening of interfacial forces in ref. [457]). The drag force influences ε G and the mean residence time of the disperse phase [420,439]. The lift force has three components, but the largest component acts perpendicular to the main bubble flow direction; in bubble columns, this means the radial direction. In particular, this force is responsible for the migration of small bubbles toward the column walls or the center of the pipe, depending on the bubble shape and size [277,278]. The turbulent dispersion force spreads out large bubbles from the pipe center and modulates peaks of small bubbles near the pipe walls [458]. Its magnitude is also high near distributor inlets, supporting the modeling of bubble dispersion near coarse gas spargers [429]. The wall lubrication force models the lift force close to the wall, thus pushing the bubbles away from it. Rzehak et al. [459] compared various formulations applied to vertical bubbly flow in a pipe and concluded that the inclusion of this force into the model is fundamental. The virtual mass force arises from the relative acceleration of an immersed moving object to its surrounding fluid; this force is often found to be negligible in bubble column simulations [421,428,437,439,449,460], and different studies have neglected it [420,430,440,441,446,450,458,[461][462][463][464][465][466][467][468][469]. Recently, Ziegenhein et al. [10] demonstrated that the virtual mass force has an influence on the prediction of the turbulence intensity at higher flow rates. It is clear that an independent validation of each single force is not possible; for this reason, a complete set of interfacial forces should be used [9,10,[443][444][445]447,459,[470][471][472][473], as discussed by [443] and recently validated against a large set of experimental data [474].
Most of the interfacial force correlations require, as an input, the average equivalent diameter of the bubbles; in this regard, two approaches can be used: (a) a constant bubble equivalent diameter (taken from experimental data or correlations); or (b) implementation of a population balance model (please, refer to the next section for further detail) that predicts the local bubble size distributions using coalescence and breakage kernels [433,[475][476][477][478][479]. The former case is typical of the mono-dispersed flow regime [480]. The latter case is typical of the pseudo-homogeneous flow regime [481,482]; in this case, a continuous approach or a discrete approach can be used. The former approach has been recently proposed by Vik et al. [292], where a continuous solution of the density function is obtained rather than a discrete solution. Moreover, both the bubble velocity, composition and temperature are all continuous functions of the bubble size. Hence, an infinite distinction is used between the small and large bubbles in an inhomogeneous population balance model. In the latter approach, the gas phase is subdivided into different bubble size classes along with a single ("homogeneous population balance model") or multiple ("inhomogeneous population balance model") gas velocity fields, depending on the desired level of distinction between small and large bubbles [465]. Homogeneous population balance models have been applied to bubble column simulations and upward bubbly flows by several authors [428,433,440,441,[461][462][463]466,475,483]. On the other hand, Krishna et al. [484] introduced one of the first uses of two velocity groups to distinguish the dynamics of small and large bubbles. This approach has been also presented in other works [483,[485][486][487]. The model, however, did not implement a population balance model, but rather, a constant bubble equivalent diameter for each group and the non-drag forces were neglected. More recently, Ziegenhein et al. [10] used two velocity groups for the gas phase, following the proposal of Krepper et al. [488], which has also been used in several studies [430,435,442,444,446,447,458,465,489]-the two velocity groups for the gas phase are computed by slitting the original distributions at the diameter for which the lift coefficient changes its sign. It is worth mentioning the study of Xu et al. [483], which compared the above-mentioned approaches and found that the best performances are given by the "inhomogeneous population balance model". Recently, Besagni et al. [481] studied the pseudo-homogeneous flow regime in large scale bubble columns and compared different numerical approaches: (i) the fixed polydispersed formulation (using the two-and four-class approaches); and (ii) coalescence and break-up closures. The results suggest that the inclusion of coalescence and break-up closures are essential to simulate the bubble column fluid dynamics.
It is worth noting that suitable boundary conditions should be provided to obtain realistic solutions of the problem-i.e., wall, inlet, outlet, and symmetry conditions. The wall boundary conditions model the interaction of the dispersed/continuous phase with a continuum material (i.e., a solid phase); at the walls, the gas shear stress and the tangential components of the liquid velocity should tend towards zero. Inlet and outlet boundary conditions aim to provide boundaries to the numerical computational domain, and they represent the effects of the neglected part of the whole domain in the considered spatial discretization. It is worth noting that, when developing a numerical model, proper selection of the boundary conditions is essential for reaching reliable numerical solutions (i.e., some reliable assumptions should be provided when modeling inlet conditions, i.e., fully developed flow conditions). It is worth noting that inlet conditions are particularly important; in practice, we generally know the superficial velocities but not the local volume fractions and the local phase velocities. Normally, the local volume fractions are estimated roughly, and the local velocities are then computed from the superficial velocity definition. It usually follows that steep gradients occur close to the inlet, as the solution of the Navier-Stokes equations in the first point predicts a solution far from the estimated profiles. In addition, it should be noted that, from a practical point of view, the convective terms in the gas phase momentum equation are generally very small compared to the pressure, drag and gravity terms. This induces convergence problems, as the gas velocity at the first point within the domain is hardly related to the boundary condition, and the finite number representation of the computer may prevent convergence (see ref. [490]). A typical outlet condition concerns a constant pressure, which requires vanishing tangential velocity components of both the gas and the liquid phases. Finally, symmetry conditions aim to truncate the numerical domain. For transient calculations, the initial conditions should be provided as well.

The Approach
An uncertainty in bubble column design and scale-up arises from the lack of understanding of the phenomena governing bubble size distributions (BSDs) and, consequently, the column fluid dynamics and the interfacial area. The bubble size distributions depend on the many phenomena occurring in the bubble column: the balance between the coalescence and break-up rate, the pressure change, the gas-liquid mass transfer. In particular, the bubble size distribution is considered, mainly as a consequence of the coalescence and break-up phenomena occurring in the bubble column. Theoretically, based on the knowledge of these phenomena, the BSD could be determined a priori. Based on the BSD of the system, the whole column fluid dynamics could be predicted. However, because of the unsatisfactory understanding of the mechanisms that govern break-up and coalescence, no broadly applicable model is available. In bubble columns, the initial bubble size distribution is determined by the formation of the bubbles at the gas sparger level. This bubble size distribution may not be stable due to changes by break-up and/or coalescence mechanisms and determines the "developed" (equilibrium) BSD existing in the bubble column. Since the interfacial area concentration changes with the variation in the bubble number density due to coalescence and break-up phenomena, a population balance model can be used to provide a statistical formulation to describe the dispersed phase in multi-phase flow [475,[491][492][493][494][495][496]. In the next sections, the fundamentals of this method and the governing equations are discussed, and then the coalescence and break-up phenomena are presented. Herein, the basis for understanding Chapter 5 (where a population balance model is formulated and implemented) is given. For further detail, the reader should refer to the review of Liao and Lucas concerning coalescence [215] and break-up [214] phenomena and modeling.

The Governing Equations
A population balance model can be described by a transport equation derived from the Boltzmann statistical transport equation. It reflects particles entering or leaving a control volume via several mechanisms. The bubble number density transport equation is often referred to as the Population Balance Equation (PBE) [215]: where n → x , V b , t is the bubble number density distribution function and specifies the probable number density of bubbles at a given time (t), in the spatial range, d( t is a source term, which reads as: where S c , S b S ph , S p , S m and S r are the bubble source/sink terms due to coalescence, break-up, phase change, pressure change, mass transfer and reaction, respectively. The coalescence source term, S c , reads: where the first term is the birthrate of bubbles of volume V b due to the coalescence of bubbles of volume, V b − V b and V b ; the second term is the death rate of bubbles of volume V b due to coalescence with other bubbles. Finally, Γ(V b , V b ) is the coalescence rate between bubbles of volume V b and V b .
The break-up source term, S b , reads: where the first term is the birthrate of bubbles of volume V b due to the break-up of bubbles with volumes larger than V b ; the second term is the death rate of bubbles of volume V b due to break-up. Finally, Ω(V b ) is the break-up rate of bubbles of volume V b , m(V b ) is the mean number of daughter bubbles produced by break-up of a parent bubble of volume V b and P(V b , V b ) is the probability density function of daughter bubbles produced upon break-up of a parent bubble with volume V b . Neglecting the interfacial mass transfer, the reaction contribution, the phase variations and the changes due to the pressure gradients, the coalescence and break-up rates are the only phenomena occurring. Considering a steady state, it follows that: Which can be written in terms of the bubble diameter: S c and S b are determined by modeling the bubble break-up and coalescence rates.

Bubble Coalescence Phenomena and Modeling
Bubble coalescence results from collision between bubbles. Attention may be restricted to binary collisions as the probability of interactions involving more than two bubbles is significantly lower in the homogeneous flow regime. Collision between two bubbles requires a non-null velocity difference between the bubbles [495]. The collision may be a consequence of at least five mechanisms ( Figure 40) [215]: (i) different rise velocity (buoyancy); (ii) mean velocity flow gradient; (iii) turbulent fluctuations of the carrier phase; (iv) bubble capture in an eddy; and (v) bubble-wake interactions. The collision between bubbles can result either in coalescence or in bouncing. Hence, to model the coalescence rate, one has to consider both the frequency of collisions (the "collision frequency") and their efficiency of causing coalescence (the "coalescence efficiency"). Coalescence models have been extensively studied over the past decades, and a variety of models have been published [215]. A first class of models is the empirical-based models, based on the modeling of the coalescence rate (without modeling separately "collision frequency" and the "coalescence efficiency") by using functions (i.e., power law) with parameters that can be adjusted to fit a set of experimental data. These models mainly depend on the experimental configuration and operating conditions of the data used to calibrate the model parameters. For this reason, the application of these models to a broader range of operating conditions is questionable. Another class of models is the "physical (or phenomenological)-based models", where the coalescence rate is considered equal to the products of the "collision frequency" and the "coalescence efficiency". The collision frequency is modeled by summing the contribution of the different collision phenomena. The "coalescence efficiency" can be modelled by three different approaches: (i) empirical approaches; (ii) physical approaches based on the "energy model"; and (iii) physical approaches based on the "film drainage model". An empirical approach is the approach velocity model of Lehr et al. [497] based on experimental data, showing that a small approach velocity leads to higher coalescence rates. This criterion is the opposite of the physical energy model of Howarth [498] which predicts coalescence when the kinetic energy of two colliding bubbles is larger than a critical value. However, the physical film drainage model, firstly proposed by Shinnar and Church [499], is the most common approach used to model the coalescence rate. The coalescence process described by this model is structured into three steps: collision between two bubbles trapping a liquid film; 2.
drainage of the liquid film; 3.
film rupture and coalescence.
The last step is reached if the contact time is larger than the drainage time (the flow fluctuations limit the bubble-bubble interaction time). Different models have been presented depending on the model hypothesis (i.e., non-deformable bubble surfaces or deformable surfaces with (i) immobile; (ii) partially mobile; or (iii) fully mobile interfaces).

Bubble Break-Up Phenomena and Modeling
A bubble moving through a continuous phase experiencesstress that acts to deform it; this stress is counteracted by the surface tension. If the deforming stress is strong enough to overcome the restoring effect of surface tension, the bubble will break-up, generating two or more daughter bubbles. A break-up model considers two sub-models: the "break-up rate" model and a "daughter bubble" distribution.
The "break-up rate" model considers the break-up mechanisms. In the case of turbulent gas-liquid systems, the deforming hydrodynamic stresses may be caused by at least four mechanisms ( Figure 41) [214]: (i) shear stresses (viscous forces); (ii) turbulent fluctuations/turbulent stresses (collision with turbulent eddies); (iii) interfacial instability; and (iv) interfacial stresses [496]. In the following text, these mechanisms are briefly outlined. Break-up due to shear stresses (i.e., when in a wake region or inside a turbulent eddy) may occur whenever the particle is subject to velocity gradients around its surface. The turbulence in the carrier phase causes pressure fluctuations alongside the interface and bubble-eddy collisions. It is worth noting that most kernels are restricted to the inertial range of turbulence, an extension to the entire spectrum of turbulence is described in refs. [493,494]. Hence, the bubble is stretched in one direction and eventually undergoes rupture. Interfacial instabilities are classified into Taylor-Rayleigh or Kelvin-Helmholtz instabilities: (i) Taylor-Rayleigh instability arises whenever a lighter fluid is accelerated into a heavier fluid; (ii) Kelvin-Helmholtz instability occurs whenever the velocity difference at the interphase of two fluids is sufficiently large. In this respect, Ishii and Kojasoy [500] provided a correlation for the maximum stable bubble size, d b max : d b max = 10 cm in air-water systems in ambient conditions (therefore, interfacial instabilities are seldom considered). Finally, interfacial stresses must be considered-as the bubble rise velocity increases, the shearing-off events gain importance, and the interfacial shear force and drag force cause the shearing-off of small bubbles at the rim of large bubbles. Likewise, when the interphase relative velocity becomes too high and thus, the interfacial force is high enough, the bubbles become unstable and start to stretch in the direction of the destroying force (i.e., in counter-current mode). Generally speaking, in most cases, the carrier phase is turbulent; therefore, only the dominant break-up mechanism due to the turbulent fluctuations is considered [214]. Break-up models should also consider the distribution function of the daughter bubble size [501]. Some examples of daughter bubble size distributions are (i) statistical distribution (i.e., normal or the uniform); (ii) empirical distribution; or (iii) phenomenological distribution [214]. In this regard, the reader may refer to the experimental high-speed imaging of Solsvik et al. [502]: they conducted in an attempt to understand the break-up cascade; in particular, they found that the assumption of binary breakup in the PBM approach is acceptable only if the initial breakup is analyzed.

Summary of the Literature Survey
In this paper, we presented a literature review on bubble column fluid dynamics at different scales and deeply analyzed the flow regimes, flow regime transitions and local and the global fluid dynamics parameters. We also presented experimental techniques for studying the global and local fluid dynamic properties, and we discussed the modeling approaches for studying the global and local bubble column fluid dynamics. Taking into account all the different results discussed, it is clear that the complete description and modeling of bubble column fluid dynamics rely on the knowledge of the phenomena at different scales. In this respect, a multi-scale approach should always be considered when proposing future experimental investigations. In the current author's opinion, the fundamental parameter for understanding and modeling the bubble column fluid dynamics is the prevailing bubble size distribution in the systems. In this regard, the bubble size distribution is imposed by the gas sparger and the operation modes (i.e., the gas flow rates) and it changes accordingly with the constraints given by the system and governs the bubble column fluid dynamics. In particular, a change in the phase properties/operating conditions would affect bubble interface properties (related to the coalescence suppression/enhancement mechanisms), which modifies the bubble size distribution and shapes, thus stabilizing/destabilizing the homogeneous flow regime and ultimately, increasing/decreasing ε G . The bubble column design parameters influence the fluid dynamics because of the boundary conditions, as well as the end-effects (i.e., wall or outlet effects). Similarly, appropriate modeling approaches should be selected by taking into account the prevailing bubble size distributions in the system.

Guidelines
Based on the literature survey proposed in this paper, some guidelines to be taken into account in future studies are proposed:

1.
Different studies have investigated, in the few last decades, either local or global fluid dynamics properties; unfortunately, the studies concerning both local and global fluid dynamics properties are still limited in number. In this respect, experimental studies should always provide a multi-scale evaluation of bubble column fluid dynamics (i.e., by studying, at least, both gas hold-up and bubble size distribution); 2.
Experimental studies should always provide detailed information concerning the main bubble column design criteria (as listed in Section 2.3.2) in order to relate the experimental results to the geometrical scale of the bubble column (i.e., "large-diameter", "small-diameter", "coarse gas sparger", fine gas sparger", etc.); in particular, information on non-dimensional bubble column diameters, aspect ratios and gas sparger openings should be always provided 3.
In comparisons involving the flow regime transition, detailed information concerning the flow regime transition criteria should be specified; in addition, authors should carefully evaluate if the flow regime transition method is suitable with the applied definition of flow regime transition; 4.
When presenting a numerical approach, sensitivity studies on (a) interfacial closures; (b) time and (c) space discretization should be always be performed; 5.
The modeling approach of the dispersed phase (i.e., mono-dispersed, bi-dispersed, PBM) should always be related to the prevailing flow regime observed in the bubble column;

Outlooks
Based on the literature survey proposed in this paper, some main comments concerning future studies are provided.
In particular, future experimental studies should be mainly devoted to the following: 1.
Proposing a precise mathematical description of the flow regime transitions, as listed and described in Section 2.2.2. This approach should take into account the role of instabilities in flow regime transitions, to support the mathematical description of the boundaries of the flow regimes; in particular, this approach should be applied to develop a comprehensive flow regime map; 2.
Performing comprehensive and multi-scale experimental studies to propose complete datasets for large-scale bubble columns (viz. gas holdup, bubble size distributions and shape data, and local flow properties); such datasets would provide a valuable basis to validate numerical approaches for scaling-up purposes; 3.
Understanding the influence of interfacial properties on the "bubble-scale" and clarifying how the "bubble-scale" influences the "reactor-scale"; 4.
Proposing a unified theory to explain all the dual effects observed in the literature (e.g., viscous liquid phases, organic compounds, inorganic compounds, etc.); 5.
Performing experimental studies concerning the influences of the operating conditions and phase properties on bubble shape in dense bubbly flow conditions; in this regard, approaches for studying the three-dimensional bubble shape should be developed; 6.
Performing a comprehensive comparison of the experimental techniques for both the global and the local fluid dynamic properties; Conversely, future numerical studies should be mainly devoted to the following:

1.
Extending the validation of the "baseline" approaches, to establish a common numerical framework for both small scale and large scale bubble columns; the validation should consider both the local and the global fluid dynamics properties in bubble columns; 2.
Extending the validation of numerical approaches to model the mass transfer phenomena in bubble columns; 3.
Extending the validation of the numerical approaches to the heterogeneous flow regime.
Author Contributions: Giorgio Besagni conducted the literature survey, post-processed the literature data and wrote the paper. Fabio Inzoli and Thomas Ziegenhein revised the paper.

Conflicts of Interest:
The authors declare no conflict of interest. Bubble volume in population balance equations [m 3 ] z i (i = 1, . . . , 5) Coefficients in the aspect ratio correlation (Equation (44)) [-]