Modeling large‐scale bioreactors with diffusion equations. Part I: Predicting axial dispersion coefficient and mixing times

Bioreactor scale‐up is complicated by dynamic interactions between mixing, reaction, mass transfer, and biological phenomena, the effects of which are usually predicted with simple correlations or case‐specific simulations. This two‐part study investigated whether axial diffusion equations could be used to calculate mixing times and to model and characterize large‐scale stirred bioreactors in a general and predictive manner without fitting the dispersion coefficient. In this first part, a resistances‐in‐series model analogous to basic heat transfer theory was developed to estimate the dispersion coefficient such that only available hydrodynamic numbers and literature data were needed in calculations. For model validation, over 800 previously published experimentally determined mixing times were predicted with the transient axial diffusion equation. The collected data covered reactor sizes up to 160 m3, single‐ and multi‐impeller configurations with diverse impeller types, aerated and non‐aerated operation in turbulent and transition flow regimes, and various mixing time quantification methods. The model performed excellently for typical multi‐impeller configurations as long as flooding conditions were avoided. Mixing times for single‐impeller and few nonstandard bioreactors were not predicted equally well. The transient diffusion equation together with the developed transfer resistance analogy proved to be a convenient and predictive model of mixing in typical large‐scale bioreactors.

aspects such as the measurement technique and feed and probe placements are not readily accounted for by a correlation.Furthermore, most mixing time correlations have been developed and fitted for singleimpeller vessels without aeration.Magelli et al. (2013) derived correlations for multi-impeller reactors, though, which are more relevant for bioreactors, but only for unaerated operation.Compartment model and CFD simulations, on the other hand, have a stronger physical basis than correlations and can incorporate the geometrical and configurationrelated details that simple correlations cannot.However, despite their power and increased accessibility with modern software and computing resources, both compartment models and CFD have some disadvantages in predicting mixing times: First, the developed model and the simulation result is case-, geometry-, and flow-field-specific.For instance, a change in the flow field due to addition or removal of impellers, change of impeller type, or strong aeration requires adjustment of the model structure and a new simulation.Second, an analytical mixing time formula would be preferable over individual simulations when deriving general results or conclusions.
Overall it would be desirable that a general model would have an understandable physical foundation like compartment model and CFD simulations do, be straightforward to use like correlations are, but would not require fitting of its parameter(s) to the validation data.An alternative to correlations and the more involved hydrodynamic models is the transient one-dimensional (1D) diffusion equation, also called the axial dispersion model (Kawase & Moo-Young, 1989;Machon & Jahoda, 2000;Pinelli & Magelli, 2000).The equation can produce a single formula for mixing time, it has an easily interpreted physical basis (axial dispersion), it can include configuration details such as feed and probe placements, and it depends only on a single parameter, the axial dispersion coefficient.
The 1D diffusion equation has received relatively little attention despite these attractive features, which is probably due to the fact that its parameter is not predicted a priori, but has been fitted instead.
The overall purpose of this two-part study was to develop diffusion equations into a general model of mixing and reaction in typical largescale stirred bioreactors.The focus of this first part is on mixing times, and the aim was here to derive a predictive formula for the sole parameter of the 1D diffusion equations, the axial dispersion coefficient, without fitting the model to the validation data.Previously developed successful 1D and 2D (two-dimensional) compartment modeling frameworks (Vasconcelos et al., 1998;Vrábel et al., 2000) were reformulated as a heat transfer resistance analogy to obtain a globally averaged axial dispersion coefficient from the impeller-wise volume flow rates, which enabled the transient 1D diffusion equation to predict mixing times.A large and diverse set of over 800 measured mixing times obtained in reactor setups relevant for biotechnology was collected from literature and used to challenge the model.The model's theoretical predictions were also studied and compared to the collected literature data, and various mixing time measurement techniques were interpreted and unified in the context of the diffusion equation.In Part II of this study (Losoi et al., 2023), the focus is on 1D steady-state diffusion equations with first-and zeroth-order kinetics, which were developed to predict and characterize the potentially heterogeneous profiles of substrate, pH, oxygen, and temperature in large-scale stirred bioreactors.

| Mixing time data
A comprehensive set of mixing time data was collected from 23 published articles representing 102 different reactor setups and 832 reported mixing times.In this context a reactor setup is considered a unique combination of reactor geometry, impeller type or placement, and working medium.Table 1 shows the number of mixing times obtained from the different studies, and also divides them according to configuration and operating conditions.The mixing times were typically a mean of three to four measurements, and the working media were mostly water or glycerol solutions of varying viscosity (solutions of different strength were treated as separate media).The data covered T A B L E 1 Distribution of mixing time data obtained from literature references.

Reference
M-u M-a M-f M-t M-p S S-p Total involved aeration in dispersed, loading, or flooding regimes.Table 2 summarizes the most relevant characteristics of the collected data, which are detailed in Supporting Information: Section S2.3 and are fully available (Supporting Information: File S1).In total 472 of the times were obtained in lab-scale (liquid volume V ≤ 0.1 L m 3 ), 201 in pilot-scale In some cases, the original publications did not explicitly provide all the details that were necessary for this work such as the gas holdup and impeller power loss due to gas flow and tight impeller spacing.In these cases, the values were either estimated directly from literature or by using published correlations.All these assumptions have been marked in Supporting Information: Section S2.3 and also in Supporting Information: File S1.

| Goodness-of-fit metrics
Two coefficients of determination based on absolute and relative error, respectively, were used to evaluate the mixing time predictions of the model.In the following, f is the predicted value and y the true experimentally determined value from literature.Mean values are denoted with an "m" subscript (e.g., y m is the mean of experimental values).To facilitate comparison of the model performance with other published works, the commonly used mean relative error and coefficient of variation (relative standard deviation) 2 m were also calculated.
The conventional coefficient of determination 2 is based on the sum of squared residuals normalized by the overall variability of the experimental data around their mean.R 2 measures goodness-of-fit in absolute terms.However, the mixing time data to be predicted are both strictly positive and have an orders-of-magnitude range from 3.20 to 1840 s, which makes a statistic based on absolute error nonoptimal.A metric based on relative error would be preferred for such data.Logarithmic error f y ln( ∕ ) is a suitable measure as the logarithm deals with relative errors such that for example, both a −50% error in f y = (1 − 0.5) and a +100% error in f y = (1 + 1) have the same magnitude (ln(2) = −ln(1∕2)) (Tofallis, 2015).For strictly positive mixing times it makes sense to penalize predictions half or double the true value equally.An analogous coefficient of determination based on squared logarithmic error was therefore used as a supplementary statistic: (1) To retain similarity with the conventional R 2 , geometric mean was used in the denominator to minimize the denominator sum just as arithmetic mean is used in R 2 to minimize its respective denominator sum (Tofallis, 2015).Like with R 2 that has a maximum value of 1 and desirable values over 0, a perfect fit . The error term in Q 2 was also decomposed to quantify the systematic and random error contributions separately (Supporting Information: Section S3).

| Estimation of model uncertainty
The main parameters of the model were obtained directly from literature correlations, which involved some uncertainty.The error expected in model prediction f purely due to the inevitable uncertainty in its N parameters x i was estimated by the first-order propagation-of-error formula assuming zero covariance between the parameters: i , where σ i is the standard deviation of i.The derivatives f x ∂ ∕∂ i were calculated numerically with a centered difference.
T A B L E 2 Overview of the 832 experimental mixing times and operating conditions obtained from literature (Table 1).

| Transient 1D diffusion equation
Mixing across the working height H (m) was modeled here with the transient 1D diffusion equation: with closed boundaries ( u z ∂ ∕∂ = 0 at both bottom z = 0 and top z H = ) (Kawase & Moo-Young, 1989;Machon & Jahoda, 2000;Pinelli & Magelli, 2000).u is the normalized concentration of the added substance or tracer (initial value 0, spatial mean 1), t time (s), d the axial dispersion coefficient representing both convective and diffusive turbulent flows (m s 2 −1 ), and z axial coordinate (m).The working height includes the gas holdup α G .Solution to Equation (2) with closed boundaries and an impulse initial condition at z 0 can be found in heat transfer textbooks (e.g., Cole et al., 2010): (3) The above equation was used here to predict both tracer curves and mixing times according to their various definitions.

| Mixing time
The sum's first time-dependent term in Equation (3) dominates the solution as equilibrium is approached, and Equation ( 3) is simplified to from which the time can be solved.More terms and a numerical solution of mixing time are required if the mixing time's heterogeneity level   u 1 − is high or if either the feed or measurement point, z 0 or z, respectively, is close to H 0.5 (Supporting Information: Section S4).In most cases the one-term Equation (4) is sufficient.Sections 3.2.1,3.2.2, and 3.2.3present how the diffusion equation applies to mixing times measured with a single probe, multiple probes, or a colorimetric method, respectively.

| Single probe
The most of the data collected in this work were obtained with a single probe measuring conductivity (586 mixing times out of 832), which increases linearly with the local concentration of a salt solution tracer.For a single probe at z, the mixing time is readily solved from Equation (4): (5) Heterogeneity level of the mixing time is specified by   u 1 − .The most common heterogeneity levels are 5% and 10%, which correspond to u = 0.95 and u = 0.90 when equilibrium is approached from below (probe far from tracer's injection point) or to u = 1.05 and u = 1.10 when from above (probe close to injection).For conve- nience, the absolute value of the logarithm's argument can be used such that u in Equation ( 5) is the homogeneity level between 0 and 1 regardless of whether the actual normalized signal (u in Equation 3) rises or decays to 1.

| Multiple probes
In some studies, multiple probes have been used, and the final mixing time can be the mean of each probe's individual mixing time (e.g., Bernauer et al., 2022;Xing et al., 2009) or the mixing time determined from an averaged signal of the probes (Mayr et al., 1992).In such cases, it is straightforward to first calculate separate mixing times using Equation 2 .For best comparison with experimental data, axial coordinates should be repeated to match the number of multiple probes occupying the same height if the probe numbers differ between axial locations.With N probes the mixing time becomes (Supporting Information: Section S5) when the feed point is not too close to H 0.5 .In numerical cases, the whole volume is usually monitored and the continuous definition is most appropriate, which yields (Supporting Information: Section S5): assuming again that the feed point is not at H 0.5 .Interestingly, a symmetric placement of discrete probes at z H 6) to Equation (7) (Supporting Information: Section S5).

| Colorimetric measurements
Mixing times are also fairly commonly measured with colorimetric methods.In the starch-iodine-thiosulphate method, the complete decolorization of the vessel contents signals the mixing time.In the diffusion equation's context, such measurements can be represented by monitoring the normalized concentration at the point furthest away from the feed (last to receive the sufficient amount of the decolorization agent).The homogeneity level is related to the stoichiometric excess of thiosulphate.For example, Cronin et al. (1994) used 25% excess of thiosulphate, which means that a u = 1∕1.25 = 80% concentration is required to completely decolorize the starch-iodine-complex in the axial point furthest away from feed ).The mixing time can then be calculated with Equation ( 5).
With an inert dye the standard deviation of local mean gray values of the experiment's video recording can be monitored (Gabelle et al., 2011).Assuming that the local mean gray value is linear with respect to local dye concentration, the mixing time is essentially a standard deviation-based mixing time (Equations 6 and 7).Quantification based on pH-indicators is discussed in Supporting Information: Section S8.

| Axial dispersion coefficient
To calculate the mixing time as described in Section 3.2, the axial dispersion coefficient d is required as a parameter.Based on classical turbulence theory, it is defined as: where U is the axial velocity fluctuation (m s −1 ) and X is the integral length-scale of turbulence (m) (Kawase & Moo-Young, 1989).
Section 3.3.1 extends Equation ( 8) to a global axial dispersion coefficient using a transfer resistance analogy, and Sections 3.3.2and 3.3.3define the velocity fluctuation and length-scale terms, respectively.The literature correlations that were used in calculating the dispersion coefficient and thus the mixing times are compiled in Table 3 along with their uncertainties.The following dispersion coefficient calculation method is used also in Part II of this study for characterization of substrate, temperature, dissolved oxygen, and pH profiles in large-scale bioreactors (Losoi et al., 2023).

| Transfer resistance analogy
Here, in this study, the structure of previously published successful and predictive 1D and 2D compartment models (Vasconcelos et al., 1998;Vrábel et al., 2000) was formalized into a resistances-in-series model analogous to basic heat transfer theory (Figure 1).The previous where H is the working height and A the tank's cross-section (m 2 ).
As illustrated in Figure 1, the transfer within each impeller stage is slowed down by a circulation resistance R C , and the transfer between impeller stages by an interstage resistance R I .
Resistances of both types are connected in series such that the total resistance is the sum of all circulation and interstage resistances: R R = ∑ i i .In this analogy, a reactor equipped with N i impellers has N i circulation resistances and N − 1 i interstage resistances.The concept of circulation and interstage resistances is coherent with the experimental findings that a smaller number of impellers in a high aspect ratio reactor results in a smaller mixing time (Cui, van der Lans, Noorman, et al., 1996;Vasconcelos et al., 1995): decreasing the number of impellers increases dispersion coefficient in the model as interstage resistances are removed.
Using heat transfer terminology, the circulation resistances are analogous to conduction resistance within solid bodies, and they are related to local dispersion coefficient d i within respective impeller regions by where H i is the considered impeller's working height.The impeller- wise dispersion coefficients are further decomposed to mechanical and pneumatic components such that d UX = i (Equation 8) is applied separately to both components.The product of cross-section A and a velocity fluctuation U is essentially a volume flow rate v, which yields where v C0 and v CG are the mechanical and pneumatic circulation flows (m s 3 −1 ), respectively, and X 0 and X G their length-scales (m).The impeller working heights are determined such that the boundaries are midway between neighboring impellers.Similarly, the interstage resistances are analogous to contact resistances between solid bodies, and they are inversely proportional to the mechanical and pneumatic exchange flow rates: The flow rates in interstage resistances were averaged from the adjacent two impellers.
Some operating conditions (impeller placement, strong aeration) do not conform to the standard model of N i circulation resistances and N − 1 i interstage resistances.Depending on the placement of the impellers, a stagnant upper zone may form (Cronin et al., 1994;Magelli et al., 2013;Vrábel et al., 1999).The size of the stagnant upper zone varies somewhat in literature, but here the working height of the top impeller was allowed to extend at most T 0.75 (Magelli et al., 2013) above the impeller itself, where T is the vessel diameter (m).The possible extra space was then considered a stagnant upper zone separated from the top impeller region by an extra interstage resistance.Furthermore, the mechanical circulation flow rate of the stagnant upper zone was set to half the top impeller's mechanical circulation flow rate (Vrábel et al., 1999(Vrábel et al., , 2000)).On the other hand, a merging flow was reported in some reactors with very tight impeller spacings.In such cases, the interstage resistances between merged impellers were removed such that only the circulation resistances within impeller stages remained (65 mixing times out of all 832).Similarly, the two bottommost impeller regions were merged and the interstage resistance between them was removed when flooding conditions were indicated in the original publications (70 mixing times out of 298 aerated).For reference, Alves and Vasconcelos (1995) extended the applicability of the Vasconcelos et al. (1995) 1D compartment model into the flooding regime by merging compartments from the two impeller regions closest to bottom.However, they also augmented the flow rates by fitting, which was not done here.

| Velocity fluctuations and volume flow rates
The mechanical circulation and interstage volume flow rates required in Equations ( 11) and ( 12), v C and v I , respectively, can be correlated to stirrer rate n (s ) −1 and impeller diameter D (m) through respective dimensionless flow numbers K C and K I : A velocity fluctuation U (Equation 8) is obtained by dividing a volume flow rate by the cross-section of the tank.
According to measurements in 0.292-0.720m tanks with Rushton turbines, the interstage flow number is (Vasconcelos et al., 1998(Vasconcelos et al., , 1995)): where the correction factor F I is approximately 1 in turbulent flow with high Reynolds numbers.P P ( ∕ ) G is the gassed-to-ungassed power ratio.Here a linear dependency on the gassed-to-ungassed power ratio was assumed in accordance with Vasconcelos et al. (1998Vasconcelos et al. ( , 1995)).The coefficient 0.20 is the mean reported by Vasconcelos et al. (1998Vasconcelos et al. ( , 1995) ) with COV≈10%.Vasconcelos et al. (1996) measured and reported the interstage flows at Reynolds F I G U R E 1 Application of the transfer resistance analogy developed here in a standard geometry vessel stirred with a Rushton turbine (radial flow) and an upwards pumping pitched blade turbine (axial flow).The overall axial dispersion coefficient d (m s 2 −1 , Equation 9) is inversely proportional to the sum of the three resistances in series (two circulation resistances and an interstage resistance).Both circulation resistances R C (s m −3 , Equation 11) are proportional to the respective impellers' working height H i (m), but inversely proportional to the mechanical (0) and pneumatic (G) circulation flows v C (m s 3 −1 ) and their respective length-scales X (m).The interstage resistance R I (s m −3 , Equation 12) is inversely proportional to the mechanical and pneumatic interstage flows v I (m s 3 −1 ).
numbers down to 200, and their results can be interpolated by setting (Supporting Information: Section S6) where nD ν Re = ∕ is the impeller Reynolds number (ν is the kinematic viscosity, m s 2 −1 ).Equation ( 16) is obviously restricted to Re > 147.
Both direct velocity measurements and 1D compartment model  15) and ( 16) were used here for all impeller types.
The circulation flow number was correlated as: by Vasconcelos et al. (1998Vasconcelos et al. ( , 1995) ) using Rushton turbines, where the correction factor F C is approximately 1 in the turbulent regime.
Equation ( 17 Section S6) which is restricted to Re>161.Similarly to the interstage flow number, Equations ( 17) and ( 18) were used here for all impeller types.The effect of gas flow on both circulation and exchange flows was based on the specific power lost and gained through aeration.
Both circulation and interstage flows were reduced in direct proportion to the impeller power loss (Vasconcelos et al., 1998(Vasconcelos et al., , 1995)).
Unfortunately, there are no direct measurements available for calculation of the gas-induced flows required in Equations ( 11) and ( 12): the induced flows reported by Vasconcelos et al. (1998Vasconcelos et al. ( , 1995) ) were fitted, not measured.Based on dimensional analysis, the ratio of pneumatic and mechanical interstage flows was assumed to be proportional to the cubic root of the ratio of pneumatic and mechanical power inputs (Vrábel et al., 1999(Vrábel et al., , 2000)).The gasinduced interstage flow was also reduced in proportion to the crosssection occupied by impellers, though linearly here for simplicity and not quartically like in Vasconcelos et al. (1998).The mechanical interstage flow rate's correction factor for low Reynolds numbers and power-reduction due to aeration were discarded.The ratio was then where gU ϵ = G G is pneumatic specific power input (W kg −1 ), g = 9.81 m s −2 gravitational acceleration, U G superficial gas velocity (m s −1 ), and ϵ L mechanical specific power input (W kg −1 ).In the absence of data, the pneumatic circulation flow was simply assumed to be equal to the pneumatic interstage flow: v v = CG IG .

| Integral length-scales
The tank diameter T has been suggested as the relevant length-scale for mixing time calculations (Nienow, 1997).The successful application of compartments with T ∕3 height (Alves et al., 1997;   Vasconcelos et al., 1995) implies that the integral length-scale would be X T = ∕3.However, this result was obtained in tanks where H T ≥ i .
Cui, van der Lans, Noorman et al. (1996); Vrábel et al. (1999Vrábel et al. ( , 2000) ) defined their predictive 2D compartment models where H T < i with three compartment rows per impeller, which indicates X H = ∕3 i instead.To accommodate both of these definitions here, their harmonic mean was used  5) and ( 7), a single-probe mixing time agrees with a standard-deviation-based one if the probe is placed at either z H = ∕4 or z H = 3 ∕4.With the probe and feed points as wide apart as possible, the most commonly used single-probe 95% and 90% mixing times (  u 1 − = 5% and   u 1 − = 10%) would be 10% and 13% higher than the corresponding standard deviation mixing times, respectively.However, if the feed is at the middle or very close to it as in Gabelle et al. (2011), the second time-dependent (k = 2) term dominates the diffusion equation's solution (Equation 3), and the two methods agree when the probe is placed at z H = ∕8 or z H = 7 ∕8 instead.In accordance with this prediction, Gabelle et al.
(2011) measured practically equal mixing times with the two techniques when the conductivity probe was located below the lower impeller, which was placed at z H = ∕6.The exact placement of the probe was not reported, but the z H = ∕8 prediction agrees well with the reported z H < ∕6 configuration.
Both Vrábel et al. (1999) and Guillard and Trägårdh (2003) reported mixing times in the same 30 m 3 reactor with an approxi- mately 22 m 3 working volume stirred with four Rushton turbines.

| Theoretical predictions
Given that the diffusion equation could coherently unify mixing times determined with differing experimental methods, it was then used to study the effects of nonideal pulses and probes on mixing times (Section 4.2.1).The resistances-in-series model (Section 3.3) was used to predict how transition from turbulent flow to lower Reynolds numbers (Section 4.2.2) and the number of impellers (Section 4.2.3)affect mixing times.

| Nonideal pulse and probe effects
According to the diffusion equation, the dimensionless mixing time becomes a rising function of the stirrer rate or Reynolds number in turbulent flow both with a finite-duration tracer pulse and a finite probe response time-constant (first-order kinetics).The whole analysis is presented in Supporting Information: Section S9, but in each case the measured dimensionless time could be expressed as In good agreement with the experimental findings, a uniform line addition ranging from the top to the middle (comparable to Figure 6b by Kasat and Pandit [2004]) results in a 12% lower mixing time, and an addition ranging one-third from the top results in a 5% lower mixing time (Supporting Information: Section S9.3).

| Effect of Reynolds number
The product of stirrer rate and mixing time, the dimensionless mixing time nt, is usually considered constant in the turbulent regime.
However, a rising trend in nt as a function of stirrer rate n has been reported in some cases at high Reynolds numbers (Gabelle et al., 2011;Guillard & Trägårdh, 2003;Rosseburg et al., 2018), and even a negative exponent a has been mentioned in the nt ~Re a functionality (Guillard & Trägårdh, 2003).An increase in dimensionless mixing time Jahoda & Machon, 1994;Vasconcelos et al., 1996).Using the flow number correction factors (Equations 16 and 18) interpolated from Vasconcelos et al. (1996), the developed model was well applicable even down to Re= 200 (Figure 2a).According to the flow-number- corrected model, approximately 2-, 4-, and 10-fold dimensionless mixing times compared to turbulent regime are found at Reynolds numbers of approximately 600, 300, and 200, respectively.i fitted for Rushton turbines (Vasconcelos et al., 1995) is shown for comparison.The data, model, and correlation predictions have been normalized by the respective values at aspect ratio 2. PBTD, pitched-blade turbine (down); PBTU, pitched-blade turbine (up); RT, Rushton turbine.Vasconcelos et al. (1995) agreed excellently with the radial flow impeller configurations to which it was originally fitted.

| Tracer curve and mixing time predictions
The diffusion equation and the resistances-in-series model for dispersion coefficient were found to be coherent with various mixing time definitions and their theoretical predictions agreed with the available experimental data.It should be noted, however, that the results in previous Sections 4.1 and 4.2 were independent of the actual values of the dispersion coefficient.The dispersion coefficient calculation procedure developed here (Section 3.3) was next validated by predicting tracer curves and the large set of experimental mixing times from literature (Table 1).Cui, van der Lans, Noorman, et al. (1996); Vrábel et al. (1999) published tracer curves measured in a large-scale reactor with and without aeration (V ≈ 22 L m 3 ).Excellent agreement was found between the experiments and the curves predicted here (Figure 3).Previous studies have shown that the diffusion equation can be fitted to tracer curves (Machon & Jahoda, 2000;Pinelli & Magelli, 2000), but in this study, the curves were predicted without parameter optimization.Sections 4.3.1 and 4.3.2discuss the multi-and single-impeller mixing time predictions, respectively, and Section 4.3.3concludes by evaluating the overall performance of the model.

| Multi-impeller mixing times
The unaerated turbulent multi-impeller data with non-pH-based measurement methods included 61 configurations (Figure 4a) and the nonflooding aerated turbulent data 20 configurations (Figure 4b).
Due to the higher variability of pH-based mixing times (Supporting Information: Section S8), they are presented and discussed separately below.Most of the data were obtained with two to four impellers (radial, axial, or a combination) in a standard geometry (H NT = L i ) with symmetrical impeller placement or close to it.The quality of the predictions was notable, and as could be expected, unaerated data were predicted better than aerated data (MRE = 18% vs.

MRE = 20%
).The few poorly predicted outliers in these data can be attributed to exotic or non-standard configurations: Some of the tight impeller spacing data with a merging flow by Magelli et al. (2013);Xie et al. (2014) were not predicted correctly by removing the interstage resistances, and the Gabelle et al. (2011) data were obtained with an unusually low impeller placement and tracer pulse exactly at the middle, which is a sensitive point in the diffusion equation's context (Supporting Information: Figure S1B).These deviant data have been annotated in Figure 4a,b.The unaerated and aerated predictions were evaluated both with and without these outliers (Table 4), and the overall performance was remarkable especially when these untypical data (28 unaerated, 10 aerated) were not considered: R ≥ 95% 2 and Q ≥ 90% 2 were achieved even in aerated data and the MRE was only 12% for the unaerated data and 18% with aeration.
The approximately normal distribution of logarithmic error also indicated a high quality of prediction (Figure 5a,b).
Flooding condition data included seven configurations (Table 1), and these data were the most poorly predicted subset by all metrics (Figure 4c, Table 4), which is also seen in the far-from-normal distribution of logarithmic error (Figure 5c).This was expected, however, for the experimental mixing times were also much less reproducible in flooding conditions (Alves & Vasconcelos, 1995;Shewale & Pandit, 2006).The merging of the two bottommost impeller regions as explained in Section 3.3.1 was insufficient to predict mixing times in flooding conditions where a bubble column like flow field starts to emerge (Alves & Vasconcelos, 1995;Shewale & Pandit, 2006).
Transition regime data were obtained in 13 configurations with Re < 10, 000 (Table 1), and they were predicted with high accuracy  and precision (Figures 4d and 5d).Some data from Cronin et al. (1994) included aeration as well, and the Alves et al. (1997) data sampled systematically various feed locations.Of all the subgroups shown in Table 4, these data were predicted with the highest R 2 and Q 2 and lowest MRE.All the data in this group were obtained with Rushton turbines, mostly D T = ∕3 in size.Two of the Magelli et al.
(2013) configurations in this group had 8 and 12 impellers in a H T = 4 geometry where four impellers would be expected, but good predictions were nevertheless obtained by the removal of the interstage resistances in the 12-impeller configuration where merging flow was reported.All the other configurations had the standard aspect ratio H NT = L i .Three studies (Table 1) reported pH-based multi-impeller data that were obtained in nine different large-scale configurations with working volumes from 1.8 up to 22 m 3 (Figure 6a).The impeller types were varied: Guillard and Trägårdh (2003) data were obtained with only Rushton turbines, Xing et al. (2009) with axial flow impellers, and Rosseburg et al. (2018) with combinations of both types.All of the Xing et al. (2009) data were aerated, and the predictions were precise (small random error) but very inaccurate with a large bias to low values.However, they obtained their mixing times in a bicarbonate buffer with addition of a strong base, and it is plausible that the pH changes were toward the equivalence point between K p a values.In such a case longer than true mixing times would be expected (Supporting Information: Section S8), which is equivalent to predictions being systematically too low.The 22 m 3 unaerated data reported by Guillard and Trägårdh (2003) were fairly well predicted with good accuracy and decent precision as well.The random error was, however, larger than what was obtained in the same configuration by Vrábel et al. (1999) with a linear mixing time determination method (Section 4.1).Rest of the Guillard and Trägårdh (2003) data were both aerated and unaerated, and the predictions deviated more from the experimental values.Three of these mixing times were obtained with two impellers in a low aspect ratio of only H T = 0.84 , and these outliers are annotated in Figure 6a.The Rosseburg et al. (2018) data (unaerated, aerated, and flooding), that were obtained by monitoring the mean gray value of a pH-indicator solution, could not be predicted with high quality.Their 95% mixing times could be interpreted as the time points where the bottom 5% of the reactor had a pH above 8.2 and the rest a pH below 8.2 (Supporting Information: Section S8).Unfortunately, it was not possible to determine in retrospect which normalized concentration of the added acid (u in Equation 5) corresponded to pH 8.2, and singleprobe 95% mixing times at z H = 0.05 were predicted as the best guess requiring least assumptions.The distribution of logarithmic error in pH-based mixing time predictions resembled a bimodal mixture of two normal distributions (not shown), one associated with the underpredicted Xing et al. (2009) and low aspect ratio Guillard and Trägårdh (2003) data and the other with the rest of the data.
Overall these pH-based data were poorly predicted, which is seen as negative R 2 and Q 2 and a high MRE in Table 4.However, difficulty in predicting was expected due to acid-base-chemistry's influence (Supporting Information: Section S8).The pH-based mixing times seem to be subject to chemistry-related case-and study-specific variations that cannot always be accounted for in modeling, which is regrettable, since pH-based measurements often are the only practical alternative to measure mixing times in large-scale reactors.
The chemical error of pH-based mixing times can be kept to a minimum, though, by (1) keeping the initial pH close to a K p a value of the buffer (tap water is a carbonic acid buffer), (2) making the pH change always toward the K p a value, (3) employing small pH changes (Supporting Information: Section S8).For example, a pH-change from K p + 0.25 a to K p − 0.25 a induces an error less than 3% to 90% mixing times.The time-constant of the pH probe's response should also be kept small compared to the measured times (Section 4.2.1).

| Single-impeller mixing times
The non-pH-based single-impeller mixing times (Table 1) were quantified with starch-iodine decolorization and conductivity techniques (Figure 6b) in 18 mostly pilot-scale ( V 0.1 < < 1 L m 3 ) configurations (sources referenced in Table 1).Impeller placements ranged from H 0.125 to H 0.5 and diameters from T 0.09 to T 0.45 .Here, the best predictions with only at most 10% errors on average were obtained for Khang and Levenspiel (1976); Langheinrich et al. (1998).
Rushton turbine configuration data which included both measurement methods.However, the predictions for Khang and Levenspiel (1976) small Rushton turbine data (D T ≤ 0.3 ) were in an average sense only 50%-75% of the true values, and all axial flow impeller data by Khang and Levenspiel (1976); Langheinrich et al. (1998) were vastly underpredicted.The overrepresentation of impellers with low power numbers (axial flow) in the underpredicted subset suggests that the impeller's specific power should not be neglected (Nienow, 1997).Overpredictions were found only in reactors with aspect ratios 2 (Cronin et al., 1994) and 3 (Vasconcelos et al., 1995).

| Model evaluation
The predictions yielded MRE = 26%, R = 92%  negligible and always smaller than due to random error in each of the considered data subsets (Supporting Information: Table S2).With a COV = 10% in circulation and interstage flow numbers and the slight uncertainty in gas holdup and power loss due to aeration (Table 3),  Vrábel et al. (1999Vrábel et al. ( , 2000) ) predicted unaerated and aerated mixing times in four large-scale reactors using 2D compartment models and reported MRE = 4% and COV ≤ 19%, respectively.
Comparable values were obtained here with MRE = 8% and COV = 18%.It seems reasonable to say that the model developed here has performed excellently and particularly with multi-impeller configurations, where also the error distributions indicated only little systematic error (Figure 5a-c, Supporting Information: Table S2).It is to be noted that the diffusion equation's only parameter was calculated with a predictive model with no fitting to the data.

| Future improvements
Standard multi-impeller configurations were predicted remarkably well regardless of whether they involved radial or axial flow impellers or a combination of them, but some non-standard and single-impeller configurations and aerated cases left room for improvement: (1) The circulation and interstage flow numbers with different impeller types were assumed here to be the same as was determined for Rushton turbines by Vasconcelos et al. (1998Vasconcelos et al. ( , 1995)).Especially the axial flow single-impeller data suggested, that the specific power input might be worth including in determining the circulation resistance.In accordance with the  2D compartment models by Cui, van der Lans, Noorman, et al. (1996); Vrábel et al. (1999Vrábel et al. ( , 2000)), an exchange flow could also have been included in the circulation resistances within impeller stages and not only in the interstage resistances.(2) Both mechanical and pneumatic circulation flow length-scales were assumed here to be limited by the tank diameter and the impeller or vessel working height.The choice to use their harmonic mean was successful, but arbitrary, and other formulations could have worked equally well or even better.In the cases of merged flow due to tight impeller spacing or impeller flooding, the lengthscales could benefit from revisiting.(3) The gas-induced flow was determined as an initial guess from specific power by dimensional analysis, which yielded a fair result unless the aeration rate was high enough to cause impeller flooding.The prediction accuracy and precision were, however, smaller than in unaerated data (Supporting Information: Table S2).The merging of the lowest impeller region into the next impeller region was insufficient to model the effects of excessive gas flow where pneumatic agitation was dominant.However, this is less of a concern for configurations with contoured-blade radial turbines at the bottom instead of flat-blade Rushtons, as they tolerate higher aeration rates without flooding.( 4 Section 3.1 presents the transient 1D diffusion equation and its solution, Section 3.2 discusses the determination of mixing time in the context of the diffusion equation, and Section 3.3 details the calculation of the required axial dispersion coefficient.

( 5 )
and to average them or to average the signals first (Equation 4) for mixing time quantification.The standard deviation of the local concentrations may also be tracked.In experimental cases a discrete definition of standard deviation is used: σ 1D compartment models and by extension the 1D axial diffusion equation here represent both the turbulent exchange and the convective axial-radial circulation flow patterns by a diffusive exchange flow, or axial dispersion.It was first recognized that the overall axial dispersion coefficient d is inversely proportional to an overall transfer resistance R (s m −3 ): fits by others suggest that the interstage flow numbers might be approximately the same also with axial flow impellers:Jahoda and Machon (1994) measured mixing times with multiple Rushton turbines and pitched blade turbines in both up-and down-pumping configurations and fitted very similar exchange flow numbers in a 1D model with N i compartments regardless of the impeller type.Vrábel et al. (2000) measured magnetically the axial velocity fluctuation away from the impellers in over 20 m 3 working volumes and found that the normalized velocity fluctuation induced by a Rushton turbine and an axial flow impeller was practically identical.Consequently, Equations ( ) agreed with their experimental measurements for D T = ∕3 and D T = ∕2 impellers.It is reasonable to assume that K C has at least the same 10% relative standard deviation as K I .Vasconcelos et al. (1996) fitted their model's circulation flow rate at Reynolds numbers down to 200, and the obtained correction factor is satisfactorily represented by setting (Supporting Information: Vasconcelos et al. (1996) remarked that at low Reynolds numbers (Re ≤ 200) the corrected circulation flows unphysically fell below the interstage flows.Interestingly, the Kolmogorov length-scale associated to the smallest turbulent eddies would have been approximately up to 3 mm in their reactor with Re = 200, which is relatively close to the scale of conductivity probe diameters.The unphysical fit might have been due to microscale mixing limitations in the lower transition region.In any case, such low Reynolds numbers are rare in bioreactor operation, and improvement of their correction factor was not attempted here.
) which favors the lower of the values.In a standard geometry, the whole working height H was used here instead of impeller-wise working heights H i , as pneumatic agitation tends to create a circulation loop encompassing the whole tank when it is the dominant form of agitation(Alves &   Vasconcelos, 1995;Machon & Jahoda, 2000;Shewale & Pandit, 2006).4| RESULTS AND DISCUSSIONWe first show that the diffusion equation accommodates to different mixing time definitions (Section 4.1).Theoretical results derived here from the diffusion equation and the developed transfer resistance analogy model regarding operating conditions, number of impellers, and nonideal tracer pulse and probe response are then presented (Section 4.2).As a final validation, a few tracer curves and the large body of experimental mixing times from literature were predicted with the model (Section 4.3).Finally, potential improvements to the model are discussed (Section 4.4).4.1 | Mixing time definitions in the context of the diffusion equationLiterature contains a couple of intriguing examples of the mixing time measurement technique's influence on mixing times.Gabelle et al. (2011) used both the common single-probe conductivity method and a dye-based method, where the standard deviation of a few local mean gray values of the tracer experiment's video recording was monitored.According to Equations ( Vrábel et al. (1999) measured 95% times with a fluorescent tracer, and their unaerated mixing times at four different stirrer rates can be summarized as a dimensionless mixing time nt = 292 ± 7 95 (mean ± sample standard deviation), which can be transformed to nt = 235 ± 6 90 by applying Supporting Information: Equation (S28) to the individual mixing times first.Guillard and Trägårdh (2003), on the other hand, reported in otherwise similar conditions nt = 221 ± 35 90 where the four dimensionless pH-based 90% times(251, 250, 200,   and 182) had a quite high COV= 16% when compared toVrábel et al. (1999) data (under 3%).Both the higher variability of the pH-based mixing times and the difference to fluorescence-based times warranted analysis (Supporting Information: Section S8): It turned out that in general pH-based mixing times deviate from "true" mixing times due to the nonlinear definition of pH and acid-base chemistry, and that the magnitude of this deviation is proportional to the magnitude of the pH change incurred by the measurement.In addition, the locations of both the initial and final pH with respect to the K p a value of the medium affect the direction and magnitude of the error.Guillard and Trägårdh (2003) mentioned that the acid pulses resulted in approximately 1 unit pH changes, which could induce even ±20% quantification errors in a tap-water like carbonic acid buffer (Supporting Information: Section S8).The nonlinearity of pH and acid-base chemistry effects alone seem sufficient to explain these discrepancies in pH-based data.Interestingly, Langheinrich et al. (1998) found that pH-based 90% mixing times matched starch-iodine-thiosulphate decolorization mixing times (25% excess stoichiometry, 80% mixing time at the point farthest away from feed) closely in one of their configurations but not in another one.Assuming that the probe was located at the impeller's height, the diffusion equation predicts exact correspondence for the two determination methods in their first configuration with H T = and z T = ∕3 bottom clearance of impeller (Supporting Information: Equation S28).Their second case with H T = 1.3 and z T = 2 ∕9 bottom clearance (assumed probe location) was less favorable: the decolorization method was reported to yield twice as long mixing times as the pH-method, but the diffusion equation predicts 19% shorter times instead.According to the equation, the methods would have agreed if the pH-probe were placed at z H ≈ 0.42 .However, that particular configuration was somewhat unusual with a very low impeller placement.
right-hand term represents "true" mixing and the second term is the error caused by nonideal probe or pulse.The error term is a rising function of the dispersion coefficient, which is directly proportional to the stirrer rate or Reynolds number in turbulent flow (Section 3.3).The error caused by the probe is approximately equal to the probe response's first-order time-constant assuming the timeconstant is at most 10% of the measured time, and the error caused by a finite pulse time is approximately 50% of the pulse's duration if the pulse lasts at most 10% of the measured time.In both cases, the error grows greater once the pulse duration or probe response time exceed 10% of the measured time.A greater degree of homogeneity is less influenced by these nonideal conditions as the true mixing term becomes more dominant, and vice versa, lesser degrees of homogeneity are more sensitive to the nonideal pulse and response times.Kasat and Pandit (2004) studied also the effect of tracer density on mixing times, and found that at greater densities the pointaddition of tracer stretched to a line addition.Based on their model fits the mixing times were on average up to 12% lower in turbulent flow with the highest tracer density.With the highest stirrer rate the mixing times were 6%-7% lower with the highest tracer density.A finite-length pulse can also be assessed with the diffusion equation: is actually expected at high Reynolds numbers as the measured time approaches the probe response and pulse duration times (Section 4.2.1).Most of the data collected in this work displayed essentially constant dimensionless times even at very high Reynolds numbers, though, and interestingly almost all of the contrasting data were obtained with nonlinear, pH-based measurement techniques.At transition flow regime the dimensionless mixing time increases noticeably at Reynolds numbers less than 10 4(Alves et al., 1997;

Figure
Figure 2b shows the model's prediction of how the number of impellers in a standard geometry (H NT = i ) with symmetrical impeller placement affects the dimensionless mixing time when probe and feed placements are kept as far apart as possible (both either at top or bottom).According to the model, increasing the number of impellers (and aspect ratio) from two to three and four results in 2.5and 4.6-fold mixing times, respectively.Experimental data by Cronin et al. (1994); Jahoda and Machon (1994); and Vasconcelos et al. (1995) were considered for comparison, and fair agreement was found: The model predictions and experimental data obtained with multiple radial impellers coincided, but with multiple axial impellers the three-and four-impeller mixing times were lower than predictions, on average 2.0 and 3.5 times the two-impeller values.The prediction for a single impeller was also too low (18% of two-impeller time, experimental references 24% and 45%).With a single impeller it is quite universally acknowledged that mixing time is related to power dissipation (Nienow, 1997), which is not included in the presented model emphasizing multi-impeller configurations.A correlation by
| 1069 Mixing time predictions in multi-impeller reactors.Mixing times determined with a pH-based method are shown separately in Figure6a.Note the logarithmic scaling of axes.Lab-, pilot-, and large-scale labels refer to liquid volumes under 0.1 m 3 , between 0.1 and 1 m 3 , and over 1 m 3 , respectively.The solid black line is the ideal x y = line, and the dashed black lines show 1.25 multiplicative error limits x y = 1.25 and x y = ∕1.25.Panels a-d show unaerated, aerated, flooding condition, and transition flow regime (Re < 10, 000) data, respectively.
The distribution of logarithmic error resembled a bimodal mixture of two normal distributions (not shown), one associated with the underpredicted axial flow impeller data and the other with the rest of the data.The single-impeller pH-based mixing times (Figure 6c) originated from Langheinrich et al. (1998) study, and they covered six configurations with a low impeller placement of T 2 ∕9, small impeller diameter D T = 2 ∕9, and working volumes from 72 L up to 8 m 3 .Apart from the three time points annotated in Figure 6c that were obtained in an untypical, very low aspect ratio of H T = 0.3 , the model performance was rather good.Rest of the pH-based data had an aspect ratio of 1 or 1.3.The three outlier points with very large relative errors resulted in a negative Q 2 (Table4) even though R = 71% 2 was decent.The distribution of logarithmic error did not resemble a normal distribution (not shown).Even with the three outliers the pH-based single-impeller predictions clearly outperformed the linear single-impeller predictions in terms of MRE (23% vs. 40%), but this is mostly due to the systematic underprediction of axial flow impeller mixing times in the linear subset.
832 mixing times encompassing all configurations, operating conditions, and measurement techniques (Table 4), which is a good score given the extent and diversity of the data.In addition to the often considered flat-blade Rushton turbines, the modeling covered other biotechnologically relevant reactor setups with only axial flow impellers such as pitched blade turbines, propellers, or hydrofoils (191 mixing times), combinations of a Rushton or other radial turbine at bottom and axial flow impellers above (132 mixing times), and other radial flow turbines with contoured blades (50 mixing times).The reduction in Q 2 due to systematic error was mostly T A B L E 4 Mixing time prediction statistics.
the model was calculated to have a 7%-10% COV depending on the number of impellers and impeller working heights.Considering that approximately at least a 7% prediction error is expected due to parameter uncertainty alone, the 12%-20% MRE obtained for nonflooding multi-impeller data with non-pH-based measurement techniques show good performance.For context: Magelli et al. (2013) reported 21% and 18% MRE for their two unaerated multiimpeller correlations that were fitted to the respective data containing 11 vessels.Here their data were predicted with a similar MRE = 17%.Vasconcelos et al. (1998) calculated unaerated and aerated mixing times for three dual Rushton turbine reactors with an ambitious MRE ≤ 5% using a 1D compartment model.Their data were predicted here with a higher MRE = 10%, which is still a fair result given that Vasconcelos et al. (1998) fitted their gas-induced flow parameter.
Cumulative distributions of logarithmic error q in multi-impeller mixing time predictions.Mixing times quantified with pH-based methods are not shown.Two normal distributions are shown for reference: both have the error distribution's variance, but one has zero mean and the other (shifted) has the error distribution's mean.The error distribution's mean and standard deviation are shown in each panel (q m and σ q , respectively).The proportion of data within a 1.25 multiplicative error ( f y ln( ∕ ) = ±0.223) is denoted by brackets.(a) Unaerated.(b) Aerated.(c) Flooding.Flooding has been indicated in the original publications.(d) Transition regime (Re < 10, 000).
Mixing time predictions for single-impeller reactors and pH-based data.Note the logarithmic scaling of axes.Lab-, pilot-, and large-scale labels refer to liquid volumes under 0.1 m 3 , between 0.1 and 1 m 3 , and over 1 m 3 , respectively.The solid black line is the ideal x y = line, and the dashed black lines show 1.25 multiplicative error limits x y = 1.25 and xy = ∕1.25.(a) Multi-impeller reactors with pH-based mixing times.(b) Single-impeller reactors with linear mixing times.(c) Single-impeller reactors with pH-based mixing times.
) The formation of a stagnant top zone or loop and the flow within such a zone could be investigated further.The same applies also for the merging flow configurations: after removing the interstage flows, the circulation flow numbers were simply assumed to remain the same as in standard geometry, which in some cases predicted mixing times correctly and in others incorrectly.It is likely that the circulation flow numbers are in reality affected by merging of the flow of adjacent impellers.5| CONCLUSIONSThe purpose of this two-part study was to develop simple 1D diffusion equations into a general model of typical large-scale stirred bioreactors, and this first part focused on predicting mixing times.A transfer resistance analogy to basic heat transfer theory was developed to calculate the diffusion equation's only parameter, the axial dispersion coefficient, from published hydrodynamic numbers, operating conditions, and reactor configuration.The proposed calculation of the dispersion coefficient was evaluated by collecting over 800 experimentally determined mixing times from literature such that diverse reactor sizes and configurations, impeller types and combinations, operating conditions, and mixing time definitions were included.Overall the model performed well, and the mixing time predictions were excellent in typical and biotechnologically relevant multi-impeller configurations even with aeration if flooding was avoided.Furthermore, the diffusion equation and the presented model for the dispersion coefficient could explain and unify different definitions of mixing time and theoretically predict general results regarding experimental conditions and reactor configuration.Thus, a simple-to-use mixing time predictor for large-scale bioreactors with a clear physical foundation was developed requiring only few literature correlations but no fitting.Part II of this study (Losoi et al., 2023) utilizes the validated dispersion coefficient to model and characterize the relevant variables in typical fed-batch operations.

rad N ax N hyb
Note: Linear refers to all mixing times that have not been measured with a pH-based technique.The * mark in unaerated and aerated refers to removing the data annotated in Figure4and mentioned in Section 4.3.1.Flooding conditions were indicated in original references.Transition flow regime data had Re < 10, 000.Multi-impeller pH-based and transition regime mixing times include aerated and flooding data as well.Symbols: N, amount of data points in (sub)group; R 2 , coefficient of determination; Q 2 , logarithmic coefficient of determination; MRE, mean relative error; N rad , amount of data points with radial impellers; N ax , amount of data points with axial impellers; N hyb , amount of data points with both radial and axial impellers.