Modeling large‐scale bioreactors with diffusion equations. Part II: Characterizing substrate, oxygen, temperature, carbon dioxide, and pH profiles

Large‐scale fermentation processes involve complex dynamic interactions between mixing, reaction, mass transfer, and the suspended biomass. Empirical correlations or case‐specific computational simulations are usually used to predict and estimate the performance of large‐scale bioreactors based on data acquired at bench scale. In this two‐part‐study, one‐dimensional axial diffusion equations were studied as a general and predictive model of large‐scale bioreactors. This second part focused on typical fed‐batch operations where substrate gradients are known to occur, and characterized the profiles of substrate, pH, oxygen, carbon dioxide, and temperature. The physically grounded steady‐state axial diffusion equations with first‐ and zeroth‐order kinetics yielded analytical solutions to the relevant variables. The results were compared with large‐scale Escherichia coli and Saccharomyces cerevisiae experiments and simulations from the literature, and good agreement was found in substrate profiles. The analytical profiles obtained for dissolved oxygen, temperature, pH, and CO 2 ${\text{CO}}_{2}$ were also consistent with the available data. Distribution functions for the substrate were defined, and efficiency factors for biomass growth and oxygen uptake rate were derived. In conclusion, this study demonstrated that axial diffusion equations can be used to model the effects of mixing and reaction on the relevant variables of typical large‐scale fed‐batch fermentations.

mean" mixing model utilizing the substrate distribution instead of its axial profile was presented (Maluta et al., 2020), showing that the level of spatial detail in hydrodynamically sophisticated simulations is not strictly necessary to correctly predict biomass yields.
At simplest, the compartment models are one-dimensional or 1D (Bisgaard et al., 2022), and they are essentially discretizations of a diffusion equation.The diffusion equation reproduces tracer curves measured in typical high aspect ratio bioreactors, as it is capable of representing both convective and turbulent diffusive forms of mixing in a 1D domain (Kawase & Moo-Young, 1989;Machon & Jahoda, 2000;Pinelli & Magelli, 2000).Coupled with suitable approximations to biologically relevant kinetics, the mathematics of diffusion could permit analytical solutions to profiles and distribution functions of the relevant variables in large-scale fed-batch processes.Such results could even be used to derive efficiency factors, which relate with a single number the performance at the large scale to a homogeneous situation (Delvigne et al., 2005).
The aim of this two-part study was to develop a general model of mixing and reaction in typical large-scale stirred fed-batch bioreactors using 1D diffusion equations.The first part derived a predictive formula for the axial dispersion coefficient.This second part focused on predicting and characterizing the profiles of substrate, pH, dissolved O 2 , temperature, and both gaseous and dissolved CO 2 using analytically soluble 1D steady-state diffusion equations with zeroth-and first-order kinetics.The cumulative distribution and probability density functions were also defined for the substrate.The modeling was compared against both experimental (Bylund et al., 1998;Larsson et al., 1996;Xu, Jahic, Blomsten, et al., 1999) and numerical (Larsson et al., 1996;Losoi et al., 2022;Pigou & Morchain, 2015) literature data concerning Escherichia coli and Saccharomyces cerevisiae fed-batch fermentations.The model solutions were also utilized to derive simple efficiency factor formulae for oxygen uptake and biomass growth rates.

| Experiments and simulations from literature
The large-scale experiments by Bylund et al. (1998), Larsson et al. (1996), and Xu, Jahic, Blomsten et al. (1999) were used as a reference for the modeling.The cited works reported glucose concentrations measured at top, middle, and bottom sections of the reactors and biomass concentrations.Dissolved oxygen tensions (DOTs) were also monitored with one probe at the middle (Larsson et al., 1996) or two probes at the middle and the bottom (Xu, Jahic, Blomsten, et al., 1999).Bylund et al. (1998) did not report the probe location, so it was assumed here to be in the middle as well.The control values for pH and temperature were provided in the referenced works.The liquid volumes were from 8 up to 22 m 3 .Table 1 lists relevant variables and quantities regarding the experiments.The axial dispersion coefficients and mixing times were calculated from operating conditions as in Part I of this study (Losoi et al., 2023).Bylund et al. (1998) reported ranges of stirrer and gas flow rates, and here the midpoint of these ranges was used as the operating condition.The mean substrate concentrations shown in Table 1 refer to time points with constant feeds and 20g L −1 biomass concentrations.Henry's constants for O 2 and CO 2 (Sander, 2023) and K p a values (Rumble, 2022) relevant for the cited experiments are listed in Table S1.Only temperature corrections were considered to Henry's constants and the K p a values, but they had only little effect, though.Taking for example ionic strength into account might have a more noticeable effect, but this was not attempted.
Altogether the substrate data from these references included 96 time points with three values each measured at the top, middle, and bottom of the reactors.Gas holdups were determined from the references from the reported liquid volumes and total dispersion heights, but for Bylund et al. (1998) experiments the holdups were estimated here using a correlation fitted for large scale (Vrábel et al., 2000).Larsson et al. (1996) and Xu, Jahic, Blomsten et al. (1999) reported oxygen transfer rate constants of k a = 180h L −1 for their setups, and by using the specific power and superficial gas velocity functionalities of the k a L correlations reviewed by Gabelle et al. (2011), it was estimated that k a L should have been approximately 70%-80% of that value in the Bylund et al. (1998) experiments.For simplicity, the same k a = 180h L −1 was used also for the Bylund et al. (1998) experiments.Some of the operating conditions for the Xu, Jahic, Blomsten et al. (1999) experiment were determined using literature based on the same large-scale reactor (Vrábel et al., 1999(Vrábel et al., , 2001)).The E. coli kinetic parameters determined by Xu, Jahic, and Enfors (1999) were used both for Bylund et al. (1998) and Xu, Jahic, Blomsten et al. (1999) experiments.
The large-scale simulations in 20m 3 liquid volumes by Larsson et al. (1996), Losoi et al. (2022), and Pigou and Morchain (2015) were used as a further reference.Larsson et al. (1996) reported glucose contours obtained with CFD simulations with standard Monod kinetics (fig.7 in Larsson et al., 1996).The cumulative distribution function (CDF) of their simulated substrate concentrations was estimated here by approximating the areas between the twodimensional (2D) concentration contour lines.Pigou and Morchain (2015) used a 2D compartment model and a metabolic model and provided heat maps and values of glucose concentration and also biomass concentrations (figs.9 and 7b in Pigou & Morchain, 2015).
Their results were considered here as CDFs and also as radially averaged 1D axial profiles.Our previously published results (Table 1 in Losoi et al., 2022) were obtained with a three-dimensional compartment model and standard Monod kinetics.The dispersion coefficient for Larsson et al. (1996) simulations was kept the same as for their experiments (Table 1).For Pigou and Morchain (2015) simulations a dispersion coefficient of d = 0.0659m s 2 −1 was calculated using the transfer resistance analogy concept presented in Part I of this study (Losoi et al., 2023) and the provided exchange, circulation, and induced flow rates (Appendix B in Pigou & Morchain, 2015).The reported 95% standard-deviation-based mixing time of 154 s yielded d = 0.106m s 2 −1 (eq.7 in Losoi et al., 2023) for the Losoi et al. (2022) simulations.
For the model used here, the volumetric (liquid-phase) substrate consumption rates r S (g L h −1 −1 ) were linearized into r k S = ,

S S
where k S is a first-order rate-pseudoconstant (h −1 ), a function of the substrate concentration S (g L −1 ).Using this definition, the rate-pseudoconstants were evaluated using mean consumption rates and mean concentrations.
Experimental mean concentrations were estimated here as weighted averages of the three measured concentrations such that the individual sampling locations were given weights according to the working height that they represented, although the weighing had only little effect on the mean.Simulated mean concentrations were either provided directly in the references or they could be calculated from the data.When comparing the model with experiments, the volumetric substrate consumption rate r S was assumed to equal the volumetric feed rate Q S (steady-state assumption), which was obtained from the references.
When considering previously published simulations, the kinetics used in the reference were analytically linearized and the first-order ratepseudoconstant k S was evaluated with the mean concentration of substrate,   S .Larsson et al. (1996) used standard Monod kinetics for substrate consumption, which yielded where q = 1.7gg h S −1 −1 is the biomass-specific maximal substrate consumption rate, X the biomass concentration (g L −1 ), and K = 0.18 gL S −1 the Monod constant.The kinetic parameters were Losoi et al. (2022).The substrate consumption rate in the Pigou and Morchain (2015) metabolic model was based on defining an equilibrium biomass growth rate and using a Pirt-form of biomass yield to calculate the anabolic demand of substrate.Their model included also the catabolic demand of substrate, which accounted for oxidative capacity and the state of the population.Here, the effects of acetate, oxygen, and population state were neglected, which ultimately simplified the first-order rate-pseudoconstant into the same Equation ( 1), but with and K = 0.05gL Notes: Larsson et al. (1996) conducted two experiments, one with the feed at the top (t) and the other with the feed at the bottom (b).Bylund et al. (1998) reported ranges of values for stirrer and gas flow rates (shown in parentheses here), and their means were used.Green's function method (Cole et al., 2010) was used to solve the steadystate 1D diffusion equations with zeroth-and first-order kinetics, or Laplace and Helmholtz equations, respectively.The method centers around integrating the considered equation's impulse response to a volumetric source term under the imposed boundary conditions.Symbolic computation software (sympy) was used for some of the derivations.

| Model statistics and uncertainty
Model fits were assessed against the experimentally determined data with the same two coefficients of determination that were used also in Part I of this study (Losoi et al., 2023).R 2 is the conventional coefficient of determination based on residuals f y − , whereas Q 2 is an analogous coefficient of determination defined with logarithmic error q f y = log( ∕ ).
The error term in both coefficients was also decomposed to systematic and random error components.Details on these metrics are given in Section 2.2 of Part I and Section S3 of Part I. Like in Part I, the error σ expected in model prediction due to the uncertainty of its N parameters x i was estimated by propagation of error with zero covariance between parameters: i , where σ i is x i 's standard deviation.In Part I, an error of σ d ∕ = 7% d was determined for the dispersion coefficient d.

| Software
The Python programming language version 3.8.5 (https://www.python.org)was used for all calculations and derivations with the packages numpy 1.19.2 (Harris et al., 2020), pandas 1.1.3(McKinney, 2010;The Pandas Development Team, 2020), scipy 1.5.2 (Virtanen et al., 2020), and sympy 1.6.2(Meurer et al., 2017).Both experimental and simulated previously published data were digitized from original figures with WebPlotDigiti- zer (Rohatgi, 2020) (5) was obtained by defining a dimensionless substrate concentration   u S S = ∕ normalized by the mean and a dimensionless axial is the ratio of the time-scales of mixing (H d ∕

2
) and substrate uptake S S S ).The square root of the time-scale ratio is analogous to the Thiele modulus used in chemical reaction engineering to characterize mass transfer in catalytic reactions.For example, M = 4 2 indicates that mixing is outperformed by reaction, as the rate of linearized substrate uptake is four times the rate of mixing.A general feel for the time-scale ratio and mixing limitations can be given by taking the longest t 95 (95% mixing time with probe and feed as wide apart as possible) as the measure of mixing rate according to eq. ( 5) in Part I of this study (Losoi et al., 2023),   S K = = 0.05 g L S −1 as the mean substrate concentration and Monod constant, both quite likely values in fed-batch operations (Bylund et al., 2000(Bylund et al., , 1998;;Castan & Enfors, 2002;Larsson et al., 1996;Xu, Jahic, Blomsten, et al., 1999), and q = 1 g g h S −1 −1 as the biomass-specific maximal uptake rate.Under these conditions Equation ( 6) is simplified to Another approach to evaluate the time-scale ratio is to utilize the steady-state approximation r Q ≈ S S , which gives the time-scale ratio as a function of substrate feed rate instead (assuming Tables 2 and 3 list example values for substrate time-scale ratio with some common 95% mixing times, biomass concentrations, and feed rates using Equations ( 7) and (8), respectively.With over 40gL −1 biomass concentrations or 16gL h −1 −1 feed rates, mixing limitations ) seem likely even in small-scale reactors with only t = 10s 95 mixing times.In large-scale reactors, where t > 200s 95 and longer mixing times are possible, mixing limitations may occur with biomass concentrations as low as X = 5 gL −1 or with feed rates as low as Equation ( 5) was solved with Green's function method (Cole et al., 2010), which allowed flexibility in defining the volumetric source Q S .Using a Dirac delta point source at x δ x x , ( − ) 0 0 , and insulated boundaries ( u x ∂ ∕∂ = 0 at both x = 0 and 1), the axial profile of dimensionless substrate concentration was found to be It should be noted that u depends on the square root of the substrate time-scale ratio, M. The CDF of (dimensionless) substrate concentration was found by first noting that a randomly chosen point in the reactor obeys the uniform distribution such that the CDF of the (dimensionless) axial coordinate is F x x ( ) = when x 0 ≤ ≤ 1. Solving for x in Equation ( 9) allowed identifying the substrate's CDF as 1 arcosh + arcosh max( , ) . .The substrate concentration's probability density function, also known as volume distribution (Morchain et al., 2014), was obtained by differentiating Equation ( 10) with respect to u (Section S1).The variance σ 2 of substrate concentration is found easiest by integrating spatially u x ( ( ) − 1) 2 with respect to x,

| Dissolved oxygen
The time-scales of mixing were (3.5m) /0.094 m s = 130s with spatially dependent zeroth-order kinetics, where the mean oxygen demand rate ODR (g L h −1 −1 ) was determined from the volumetric substrate feed rate using a constant yield coefficient (Bylund et al., 2000;Xu, Jahic, & Enfors, 1999): The local consumption was considered to have the same axial profile as the substrate concentration.The corresponding mass balance was where k a ), and u the local dimensionless concentration of substrate (Equation 9).
Solving for liquid-phase O 2 and limiting the values from below to zero yielded Symbols: M 2 , substrate time-scale ratio (Equation 7); X, biomass concentration; t 95 , 95% mixing time with widest possible feed-probe distance.
T A B L E 3 Example values for substrate time-scale ratio with feed-based calculation .In the referenced studies the flow of air into the bioreactors was so high that even with a 1 gg −1 consumption of oxygen per substrate the overall gas-phase oxygen conversions could have been 21%-48% at most.The effect of gas-phase depletion was then neglected, which simplified the treatment.
The simple zeroth-order formulation allowed approximating the oxygen-limited volume fraction of the reactor directly as the volume fraction where the substrate concentration induced a demand u ODR * exceeding the maximal transfer rate where the threshold substrate concentration is Since the maximum demand is found at the feed point x 0 , the gas-phase concentration O G in Equation ( 16) was evaluated at the feed point's hydrostatic pressure as well.

| Efficiency factors
Using the distribution and density functions derived for substrate, it was possible to derive efficiency factors for oxygen uptake and biomass growth on the main substrate.The efficiency factors represent the fraction of oxygen demand that the transfer rate is capable of satisfying and the fraction of substrate uptake that the population is adapted to continue growing on.At substrate concentrations below the threshold u* (Equation 16), where the local substrate-induced demand equals the transfer rate, the uptake rate equals the local demand such that u u OUR( ) = ODR , but at concentrations above the threshold the uptake rate is limited to the transfer rate such that The efficiency for oxygen uptake rate was obtained by integrating the substrate-dependent oxygen uptake rate u OUR( ) with respect to the substrate concentration and dividing this overall volumetric oxygen uptake rate by the overall where f is the substrate's density function (Equations S1 and S2), where G is the first moment of the substrate distribution (Equation S7).
The efficiency of growth on the main substrate was determined similarly by utilizing the population balance and adaptation concepts (Morchain & Fonade, 2009;Morchain et al., 2013): the population was assumed to grow at the growth rate μ u ( ) allowed by the environment where substrate concentration was below the mean (u < 1), and to grow at the mean growth rate   μ where substrate concentration exceeded the mean.Here the growth rate allowed by the environment was assumed to follow the substrate profile just like the oxygen consumption, To give an interpretation of the efficiency factors, they can both be transformed into simplistic biomass yield efficiencies.Given an aerobic biomass yield on glucose Y = 51% 0 and an anaerobic yield Y = 15% 1 (Xu, Jahic, & Enfors, 1999), the oxygen-uptake-based yield efficiency was The concept is that globally the microorganism utilizes glucose aerobically as far as possible and that the rest of glucose consumption is anaerobic.For biomass growth a similar yield efficiency can be formulated by using η μ , aerobic yield on glucose , and aerobic yield on acetate that has resulted from glucose overflow Y = 0.667 × 0.4 = 26.7% 1 (Xu, Jahic, & Enfors, 1999).For growth rate efficiency the concept is that the population grows on glucose, but dissimilates glucose to acetate by overflow metabolism when substrate concentration exceeds the mean and then consumes the acetate later on.

| Temperature
The steady-state temperature profile was determined by using a heat source with the substrate's spatial distribution and a uniform cooling that balanced the heat source across the whole volume.In effect the diffusion equation for temperature had then two zeroth-order kinetic terms, one spatially variable and the other spatially uniform.The balance for temperature was where ρ = 1000kg m −3 is the fermentation broth density, C = 4180Jkg k p −1 −1 the specific heat capacity of water (Rumble, 2022), T temperature (K), and (Doran, 2013).Equation ( 22) was solved with insulated boundaries and simplified into by defining a nondimensional temperature The substrate time-scale ratio M 2 and feed point x 0 define the shape of the dimensionless temperature profile.

| Gaseous and dissolved carbon dioxide
The profiles of CO 2 (aq) and CO 2 (g) were estimated for complete liquid-phase saturation, where all produced CO 2 is released as gas.
Plug flow was assumed for the gas phase (Dahod, 1993;Royce & Thornhill, 1991), and the mass transfer rate constant for CO 2 was considered to be 89% of O 2 's k a L (Royce & Thornhill, 1991).For the experimental references considered here, assuming the same k a L implies that the liquid-phase mixing of CO 2 (aq) is surpassed by local gas-liquid transfer just as with O 2 (Section 3.2).The overall CO 2 production was related to the substrate feed rate by Q with a constant yield coefficient and the local CO 2 production rate to the substrate profile u.The CO 2 and O 2 yield coefficients correspond to a respiratory quotient of 1.The resulting balance equation was where n D is the molar flow of CO 2 (g) (molh −1 ) and Q D the CO 2 production rate in the whole working volume divided by the crosssection (molm h −1 −1 ).Inlet molar flow at the bottom was considered to be zero, but a nonzero value could be used as well according to the composition of the fed gas.The molar flow of CO 2 (g) integrated to Assuming then that the total molar flow of gas remains constant throughout the reactor, the partial pressure and concentration of CO 2 (g), p D , and D G , can be obtained using the ideal gas law accounting for hydrostatic pressure.The local steady-state balance for CO 2 (aq) under negligible liquid-phase mixing is  (Royce & Thornhill, 1991).In accordance with experimental data (Dahod, 1993), Equation ( 27) implies that a higher stirrer rate and thus a higher k a L would decrease the CO 2 (aq) excess.The concentration of CO 2 (aq) can also be expressed as the partial pressure of CO 2 (g) that would be in equilibrium with CO 2 (aq): . After determining D L , the concentrations of HCO 3 − (aq) and CO 3 2− (aq) can be obtained from the respective acid-base equilibria.

| pH
The referenced experimental media contained 10-100 mmolL − 1 concentrations of H PO

| Substrate profiles and distributions
Instantaneous spatial profiles of substrate concentration produced by the model were compared against simulated profiles obtained with a complex metabolic model in a 2D compartment model (Pigou & Morchain, 2015) and also against S. cerevisiae fed-batch cultivations in an over 20m 3 working volume (Larsson et al., 1996).Using the published parameters of the metabolic and hydrodynamic models, the substrate time-scale ratios M 2 corresponding to the considered three profiles by Pigou and Morchain (2015) were estimated to be 2.1, 13.5, and 73.1 at 7, 9, and 15 h process times with 2.4, 5.0, and ) with the radially averaged profiles obtained by the much more complex modeling (Figures 1a and S2A).
As could be anticipated for the linearized kinetics (Figure S1), the two profiles that had a mean concentration approximately equal to and lower than the Monod constant at 9 and 15 h were associated with the higher 94% R 2 values.The volumetric feed rates and experimental substrate concentrations reported by Larsson et al. (1996)  The CDFs of substrate concentration were calculated for the Pigou and Morchain (2015) data with the same substrate time-scale ratios as earlier (2.1, 13.5, and 73.1) but without radial averaging, and also for the CFD-simulation results by Larsson et al. (1996), for which the substrate time-scale ratio was estimated to range from 9.0 to 19.2 using their kinetic parameters and operating conditions.The model produced almost identical distribution functions to Pigou and Morchain (2015) data (Figure 2a) even though the R 2 varied from approximately 0% (15 h) up to 91% (7 h), but the Larsson et al. (1996) simulation data had higher variability than the model predicted (Figure 2b), resulting in R 2 values of 43% for the top-fed and 68% for the bottom-fed case.Given the inevitable inaccuracy in estimating the CDF of the CFD-simulated substrate contour curves from Larsson et al. (1996), the model performed reasonably (Figure 2b).
The most notable discrepancies between the model here and the referenced simulations were found at the upper end of the substrate distributions, which was also expected due to the higher number of spatial dimensions in the cited works.For example, the Pigou and Morchain (2015) data at 9 and 15 h with radial averaging were both reconstructed by the model with R = 94% 2 (Figure 1a), but only were obtained by the 1D model for the CDFs of the same 2D data without radial averaging (Figure 2a).
The time evolutions of local substrate concentrations were calculated with the model for the large-scale S. cerevisiae and E. coli fermentations with up to 40 h process times (Bylund et al., 1998;Larsson et al., 1996;Xu, Jahic, Blomsten, et al., 1999).Error estimates (Section 2.4) were also calculated for the model fits.et al., 1996) with top and bottom feeds.The batch with top feeding (Figure 3a) was fitted better than the bottom-fed batch (Figure 3b), where the measured heterogeneity was less than the model suggested here.As shown in Figure 3c,d, respectively, the largescale E. coli fermentations were also well fitted by the model, but the bottom-fed 8m 3 batch (Bylund et al., 1998) better than the 20m 3 batch (Xu, Jahic, Blomsten, et al., 1999), where the measurements showed less heterogeneity than the model.As shown in Figure S5,   Larsson et al. (1996).The glucose was fed at the top.(b) Otherwise the same setup as in panel A but with the glucose feed close to the bottom sampling port (Larsson et al., 1996).(c) An 8m 3 Escherichia coli fermentation reported by Bylund et al. (1998).The glucose was fed close to the bottom port.(d) A 20m 3 E. coli fermentation reported by Xu, Jahic, Blomsten et al. (1999).The glucose was fed at the top.
unexplained (38%) was solely due to random error in this metric based on absolute error.The distribution of residuals was rather symmetrical but sharper than normal, however (Figure S6B).
The two coefficients of determination were both above 50% for the substrate time series data (Figure 3), which would be relatively low for a fitted correlation, but indicates good performance here as the model was not optimized to the data.Most of the uncertainty in model predictions was caused by the fact that the mean substrate concentration was not directly available but had to be estimated using only three experimental concentration values measured at the top, middle, and bottom of the reactors.The estimated mean values were on average 22mgL −1 with a 21mgL −1 standard deviation, which is reasonable for fed-batch operations, though.The uncertainty caused by the dispersion coefficient was small in comparison.The experimental values closest to the feed point were also rather variable, which was especially apparent with top feeds.Nevertheless, the error distributions showed good quality of fit in both absolute and relative error scales (Figure S6).Most of the error was random and not systematic in nature, lack of precision instead of lack of accuracy.

| Variance, substrate time-scale ratio, and feed points
The substrate's spatial variance was calculated as a function of the substrate time-scale ratio with Equation ( 11) and also with experimental and numerical references (Larsson et al., 1996;Losoi et al., 2022;Pigou & Morchain, 2015) for comparison.According to the model, the variance starts eventually to grow linearly with respect to the square root of the substrate time-scale ratio, as the rate of reaction exceeds the rate of mixing, which was somewhat apparent also in the experimental (Larsson et al., 1996) and numerical (Pigou & Morchain, 2015) references (Figure 4a). Figure 5b shows the estimated DOT profiles, and the proportion of O 2 -limited zones (DOT= 0 with the zeroth-order kinetic approximation) was calculated to be 22%-46% (Equation 15).The Figure 5c shows temperature profiles estimated by the model for the experimental references (Table 1) corresponding to the same time instants with 20gL −1 biomass concentrations as in Figure 5a,b with substrate and DOT profile estimates.The oxygen uptake efficiencies were estimated to be 57%-85% in these situations, leading to 0.83-1.54gL h −1 −1 oxygen uptake rates.Consequently, the gas-phase O 2 conversions would have been 6%-12%.On the basis of the modeling here, the axial profile of temperature should have been virtually homogeneous in the referenced large-scale experiments.Even with a higher 1gg −1 consumption of oxygen per substrate the estimated temperature differences would have been only up to 0.16 ∘ C, which is still negligible.For a general assessment of whether axial temperature differences could be expected to occur, temperature differences were evaluated with substrate feed rate and 95% mixing time as parameters by assuming   S x = 0.05gL , = 1 −1 0 , and According to these estimates a notable temperature difference of 1 ∘ C scale could occur at large scale (t ≥ 200s 95 ) with higher than 10gL h −1 −1 substrate feed rates, provided that oxygen transfer is not limiting (Table 4).Again, a higher than 0.446g g −1 consumption of oxygen per substrate would increase the difference accordingly.No reports of axial temperature differences were found in experimental literature, and no simulation works were found either, which implies that either there is no reason to expect any major axial differences to exist or they have been neglected.The estimations for axial temperature differences in the referenced experiments were very small, which strongly suggests the former.Unless higher than 10gL h −1 −1 substrate feed rates are utilized with sufficient oxygen transfer (Table 4), the assumption of axially constant temperature remains applicable.Local temperature differences in the proximity of the heat transfer surfaces are more likely to exist.
The profile of CO 2 (g) partial pressure was different from the others in that it was modeled by plug flow without any dispersion. (a) F I G U R E 5 Model estimations of axial distributions in the top-(t) and bottom-fed (b, in gray) large-scale cultivation experiments reported by Bylund et al. (1998), Larsson et al. (1996), andXu, Jahic, Blomsten et al. (1999).The model estimations apply to time points with 20gL −1 biomass concentration (Table 1).(a) Substrate, (b) oxygen, (c) temperature, and (d) partial pressure of CO 2 (g).
Consequently, the partial pressure in the gas phase was always zero at the bottom (Figure 5d).According to the gaseous and dissolved CO 2 profiles estimated by the model for the experimental references (Bylund et al., 1998;Larsson et al., 1996;Xu, Jahic, Blomsten, et al., 1999), the mean partial pressure exerted by CO 2 (aq) may have been relatively low, at most 35 mbar, when top feeds were utilized (Figure 6a).Bottom feeds increased CO 2 (g) and CO 2 (aq) earlier on, and over 60 mbar mean CO 2 (aq) pressure was estimated in the Bylund et al. (1998) bottom-fed case.It was also noticed that the bottom-fed batches were modeled to have a relatively homogeneous profile of CO 2 (aq) (Figure 6a).The CO 2 profiles were plausible: CO 2 (aq) concentrations at the top of the reactors were 3%-87% in excess of equilibrium with CO 2 (g) (Figure 6a), which is consistent with the experimental findings of 10%-90% excesses depending on the operating conditions (Dahod, 1993).Baez et al. (2009) measured 110 mbar dissolved CO 2 in a 5-L reactor with a 60gL −1 E. coli concentration.It was suggested by Baez et al. (2009) that the CO 2 pressure might increase up to 300 mbar at the bottom of a large reactor due to hydrostatic pressure.The performed modeling suggests that in a fed-batch process such high values are obtainable at the bottom only if the feed is at the bottom as well.This is due to the typically low inlet CO 2 content of air and the relatively little gas- phase dispersion expected in large-scale reactors (Dahod, 1993;Royce & Thornhill, 1991).
The axial pH differences due to heterogeneity in CO 2 (aq) in the experimental references were estimated to range from 0.08 to 0.31 (Figure 6b).The estimated profiles reflect how the phosphate buffers adjusted to the heterogeneous CO 2 (aq) concentrations.
Bottom-feeding led to more homogeneous CO 2 (aq) and thus pH profiles as well.The lower buffer concentrations (Table 1) used by Larsson et al. (1996) and the greater distance from the buffer's K p a resulted in larger axial pH differences.Only single pH probes were utilized in the referenced experiments, and as such direct comparison or verification of the axial pH differences could not be made.It needs to be noted that these estimates were sensitive to the plug-flow assumption and the local steady-state approximation made for CO 2 (aq) in Equation ( 27).Both the gas-phase dispersion of CO 2 (g) and liquid-phase mixing of CO 2 (aq) would smoothen and homogenize the CO 2 and consequently the pH profiles to some extent.
T A B L E 4 Axial temperature differences (°C) between top and bottom evaluated with different substrate feed rates Q S and 95% mixing times assuming   S x = 0.05gL , = 1 −1 0 , and reported an 85.5% biomass in their bottom-fed large-scale E. coli batch in comparison with a small-scale batch, and here the estimated time-averaged yield effectivities were 82% based on estimated O 2 uptake efficiency of 74% and 84% based on adaptation efficiency of 67%.Larsson et al. (1996) did not report comparisons of lab-and large-scale yields.
Even though the efficiency factors were rather simplistic, they were consistent with the yields reported in the literature.Estimation of the time-averaged η OUR for Bylund et al. (1998) experiments was rather uncertain, though.Both the gas flow rate and stirrer rate were adjusted in their experiments, but here just the middle point of the range was assumed to hold for the entire duration.For their data the k a L was also only roughly approximated from the other referenced experiments.Also the time-averaged OUR-effectivities were esti- mated using a time-independent coefficient of O 2 consumption per substrate, even though the consumption is likely to increase during a fed-batch process (Bylund et al., 2000).Furthermore, the yield losses estimated here were simplified also in the sense that the effect of maintenance was not considered.Interestingly, Maluta et al. (2020) correlated yield loss to substrate variance.Here, the simple yield losses were also related to variance through both the population balance concept but also through O 2 uptake efficiency.The relation was not linear, though, unlike in the referenced modeling work.

| Assumptions, limitations, and applicability of the model
First, it should be noted that the model was not optimized to any of the experimental or numerical references, but the dispersion coefficients were calculated directly from operating conditions using the methodology developed in Part I of this study (Losoi et al., 2023), and the kinetic parameters and mean concentrations were obtained or calculated from the referenced works.The model involved two main assumptions on the kinetics: (1) The substrate consumption was considered to be linear or first-order with respect to substrate concentration.
(2) The rest of the consumption or production rates were modeled with zeroth-order kinetics by first estimating the overall volumetric rate using a global balance and then using the substrate profile to transform this total consumption or production rate to a local rate.Hydrodynamically the major assumptions were to assume (1) turbulent axial dispersion for the liquid phase and (2) plug flow for the gas phase.The solid phase (biomass) was not distinguished from the liquid phase here.
The  cultivations, where the cell population has adapted to its environment (Morchain & Fonade, 2009;Morchain et al., 2013).Scale-down experiments have shown that the Monod kinetics do not apply at dynamic, heterogeneous conditions (Xu, Jahic, Blomsten, et al., 1999): substrate uptake rates exceeding the "maximal rate" parameter have been found, when a culture is suddenly exposed to a higher substrate concentration than what it has adapted to.The linearized kinetics alleviate this to some extent by allowing the uptake rate to exceed the conventional maximal uptake rate set by the standard Monod kinetics.The multiple uptake systems existing in the considered microorganisms were all represented here by a single rate expression, which is a conventional choice despite neglecting an inherent biological feature (Morchain et al., 2021).The net uptake rate of the various substrate uptake systems with different affinities in an adapted population is often quite well represented by a single Monod-or Michaelis-Menten-type expression both in E. coli (Quedeville et al., 2018) and S. cerevisiae (Youk & van Oudenaarden, 2009), though.
Using the substrate profile to estimate the local consumption or production rates from a global volumetric rate was a convenient choice that allowed analytical profiles to be formed for dissolved O 2 , temperature, and CO 2 .This implied a one-way coupling between substrate and O 2 and substrate and CO 2 : substrate induced consumption of O 2 and release of CO 2 , but their availability or presence did not influence substrate consumption.A two-way coupling would be more appropriate, as O 2 limitation tends to increase substrate consumption (Xu, Jahic, Blomsten, et al., 1999), but in a fed-batch context the overall substrate consumption rate eventually equals the volumetric feed rate, making the straightforward one-way coupling more applicable.In a batch reactor the twoway coupling would be more critical.It needs to be noted, however, that the approach assumed the biomass concentration to be high enough for the substrate feed rate rather than the metabolic capacity to be the limiting factor.The k a L and yield coefficients are also required as a parameter by the model and need to be estimated by measurements or correlations.The obtained O 2 profiles were more indicative of potential O 2 limitation zones that could be defined here as environments with O = 0 L , but probably less applicable as definitive profiles.There is always some upper limit to biological oxidative capacity (Szenk et al., 2017), which was not considered in this modeling work, however.It would be possible to estimate the spatial profile of dissolved O 2 using also a biological limit to the local uptake rate.This would hardly bring substantial value, when the objective is simply to detect potential O 2 limitation around feed points in a fed batch.In a batch setting the biomass-specific rates and limits to them play a more important role, when the oxygen consumption is not as limited by substrate availability.Likewise the profiles of temperature and CO 2 were relatively easy to obtain using the substrate profile.For simplicity, local limitations in O 2 transfer were not considered in determining the temperature profile.Instead, the local limitations were incorporated through the O 2 uptake efficiency factor η OUR , which affected the overall O 2 consumption and heat release rate.Using a nonuniform cooling in obtaining the temperature would have affected the profile, but most likely not in a significant amount.
As for the hydrodynamic assumptions, the use of turbulent dispersion for the liquid phase was treated and validated in Part I of this study (Losoi et al., 2023).Plug-flow assumption for the gas phase was necessary here to obtain an estimate of CO 2 profiles.The assumption could be considered reasonable in the high aspect ratio reactors studied here (Dahod, 1993;Royce & Thornhill, 1991).With lower aspect ratios it is probable that axial dispersion cannot be neglected.The gaseous and dissolved CO 2 balances also implied that the liquid phase is locally saturated with CO 2 such that all produced CO 2 is released as gas and none dissolves (Royce & Thornhill, 1991).
A high O 2 flow rate was necessary for assuming undepleted gas phase.With lower flow rates relative to theoretical maximum consumption the gas-phase conversion should be taken into account.For preliminary analyses it is perhaps easiest to use zero conversion and to note that if limitations are predicted, they are likely to occur as well.In the modeling performed here, the assumption of negligible gas-phase conversion of O 2 was not too bold (conversions 6%-12%).In some other contexts, using zero conversion throughout would be unreasonable.If the nonzero conversion cannot be assumed, the profile of gas-phase oxygen could be roughly estimated similarly to CO 2 (g) by a plug-flow equation but with spatially dependent consumption.The axial and radial differences in k a L due to impeller vicinity (Oosterhuis & Kossen, 1984) were neglected for simplicity.In a heterogeneous fed-batch setting such as here the spatial variations in k a L are not necessarily as important as in a homogeneous batch setting, where the O 2 demand by substrate is more-or-less uniform across the whole reactor.It would be possible to use a spatially heterogeneous k a L when estimating the profiles, if more precise data were available, or by correlations (Oosterhuis & Kossen, 1984).Also, the space below the sparger at the bottom can usually be expected to be poorly oxygenated (Oosterhuis & Kossen, 1984), which was not accounted for here.2 and 3).

| Implications
Heterogeneous substrate concentration profiles localize O 2 demand as well, leading to anoxic zones.Similarly the use of cofeeding strategies, where an additional high-energy substrate is supplied in low concentrations (Park et al., 2019), could be compromized by the high substrate concentrations found near the feeding points.If the substrate feed rates were more intensive, for example, 10gL h −1 −1 , and O 2 transfer were not limiting, measurable axial temperature differ- ences might be expected especially if the substrate is completely oxidized for energy.
On the basis of the various negative effects of the characterized heterogeneity, it is suggested that large-scale reactors should be homogenized more effectively.Of all the alternatives, the use of multiple feed points or at least positioning the feed in the middle instead of at the top (Losoi et al., 2022) could be the easiest to implement.Symmetrical feed placement divides the effective working height by the number of feed points N such that M N 2 −2 holds, leading to a quick decrease in the time-scale of mixing and thus heterogeneity.From the dissolved CO 2 and pH perspectives a bottom feed would be most desirable.Homogeneity in the gas phase is not achievable by substrate feed arrangements, but linear gasphase composition profiles would be found in the case of a homogeneous substrate concentration.

| CONCLUSIONS
The aim of this two-part study was to comprehensively model largescale stirred bioreactors using 1D diffusion equations.Part I of this study (Losoi et al., 2023) presented a computation formula for the model's parameter, the axial dispersion coefficient, and validated it against a large set of previously published experimental data.This second part employed the model to characterize substrate, pH, oxygen, CO 2 , and temperature profiles in typical fed-batch contexts.
The characterizations were compared with available experimental and numerical data, and good accordance was found even though the model was not optimized to the reference data.The modeling suggested that indeed each of the five variables could be heterogeneous, though temperature not as severely as the others.
According to the model, appropriate feed point placement could effectively homogenize the liquid phase.CO 2 could not be homogenized in a tall reactor, but a linear profile of gas-phase content could be expected if the substrate were homogeneously distributed.Likewise, gas-phase O 2 conversion would be expected to be linear in a tall but homogeneous reactor.Using a bottom feed could improve O 2 consumption and yield more homogeneous profiles of dissolved CO 2 and pH.On the basis of this two-part study, 1D diffusion equations can be applied for simple and predictive preliminary modeling of typical large-scale stirred bioreactors.

AUTHOR CONTRIBUTIONS
Pauli Losoi developed the model, performed the computations and analysis, and wrote the manuscript.Jukka Konttinen and Ville Santala supervised the study and revised the manuscript.

ACKNOWLEDGMENTS
Financial support by Tampere University of Technology Graduate School is acknowledged.The work presented in this article is supported by Novo Nordisk Foundation grant NNF22OC0079579.
Figure S1 compares the analytical substrate profile with linearized kinetics to profiles determined numerically by finite-volume discretization and with standard Monod kinetics assuming the Monod constant is 0.1, 1, or 10 times the mean substrate concentration.The analytical profile is remarkably close to the numerically solved unsimplified profiles, yielding mostly deviated substan- tially from the linearized analytical profile, but this case corresponded to an unlikely situation having both a considerable mixing limitation (time-scale of mixing 16 times the time-scale of reaction) and a high mean concentration of substrate.Thus, in the context of the steadystate 1D diffusion equation, linearization is a good approximation to Monod kinetics provided that the mean concentration of substrate is known or can be predicted.
et al. (1998) and Xu, Jahic, Blomsten et al. (1999) data, respectively, whereas the k a = 180h L −1 corresponded to a 20 s transfer time-scale surpassing the liquid-phase mixing.The profile of dissolved O 2 was then obtained by a steady-state approximation without mixing assuming that local transfer and consumption rates are equal.The local consumption rate was estimated 27) where h D is Henry's constant for CO 2 (mol mol L G −1 ) and 0.89 the coefficient used to relate the k a L values of CO 2 and O 2 .Equation (27) allows estimating the concentration of CO 2 (aq) marked by D L (mol m −3 ), which is in excess of the local gas-liquid equilibrium h D D G to balance the local transfer rate with the local production rate

4
| RESULTS AND DISCUSSION Substrate profiles and distributions produced by the model are first presented and discussed in Section 4.1 along with both experimental and numerical reference data from literature.The effect of substrate time-scale ratio and feed point number and placement on the substrate's volumetric variance were then calculated (Section 4.2).Instantaneous profiles of O 2 , temperature, CO 2 , and pH correspond- ing to 20gL −1 biomass concentrations were estimated for the referenced large-scale experiments (Section 4.3).Finally, biomass yield effectivities were determined for the experimental references by first calculating oxygen uptake and adaptation efficiency factors directly from the experimental substrate concentration data (Section 4.4).The assumptions, limitations, and applicability of the model are evaluated in Section 4.5, and implications of the characterized profiles and distributions are discussed in Section 4.6.Section S2 contains supplementary general-level results and discussion.
-scale of mixing being twice and over 10-fold the timescale of reaction.The profiles yielded by the model were in

.
yielded substrate time-scale ratios ranging from 17 to 30.The biomass concentrations in these data were within 10-20g L −1 .The top-fed profile with 10gL −1 biomass concentration was well represented with R = 94% 2 The other three model profiles were similar despite having an unsatisfying R < 0 2 , but this could for the most part be due to comparing against only three reported experimental concentrations and the uncertainty in estimating the mean concentration from only three values.

F
the data and model estimates covered a wide range of substrate concentrations relative to the Monod constants K S .Altogether, I G U R E 1 Axial profile of substrate concentration.(a) The marks show data obtained by radial averaging from a numerical two-dimensional study byPigou and Morchain (2015).The mean concentrations were 274, 58.1, and 13.6 mg L −1 at 7, 9, and 15 h, respectively.(b) Experimental data byLarsson et al. (1996).Labelst1, t2, b1, and b2 refer to top-(t) and bottom-fed (b) runs with biomass concentrations of 10, 15, 10, and 20gL −1 , respectively.The corresponding calculated mean concentrations were 17. 7, 17.2, 16.2, and 10.2mg L −1 .approximately one-half of the model values were within one estimated modeling error from the measured value (131/288).The logarithmic coefficient of determination was Q = 52% 2 with 2% and 46% contributions by systematic and random error, respectively, to the fraction of variance unexplained.The distribution of logarithmic error was approximately normal (Figure S6A).The conventional coefficient of determination was slightly higher at R = 62% 2 , and the systematic error was negligible such that the fraction of variance (b) (a) F I G U R E 2 Substrate concentration's cumulative distribution functions.(a) The marks show data obtained from a numerical two-dimensional study byPigou and Morchain (2015).(b) The marks represent data obtained from computational fluid dynamics simulations byLarsson et al. (1996) with top (t) and bottom feeds (b) with biomass concentrations 15 and 20gL −1 , respectively.
Model predictions and measured glucose concentrations in large-scale fed-batch fermentations.The gray regions denote the estimated error of model prediction.(a) A 20m 3 Saccharomyces cerevisiae fermentation reported by Figures 5 and 6 show profiles of DOT, temperature, gaseous and dissolved CO 2 , and pH estimated by the model for the experimental references at time points corresponding to 20gL −1 biomass concen- trations and a constant feed (parameters in Table1).The calculated substrate time-scale ratios based on the estimated mean concentrations and reported volumetric feed rates were between 16.8 and 30.2, and the corresponding substrate profiles are shown in Figure5afor reference.

F
I G U R E 4 Volumetric variance of dimensionless substrate concentration.The square root scaling of the horizontal axis emphasizes the eventually linear σ M 2 relation.(a) Experimental data by Larsson et al. (1996) (the same data sets as in Figure 3a,b).Simulation data by Pigou and Morchain (2015).(b) Simulation data by Losoi et al. (2022).bottom-fed cultivations showed less limitation due to higher hydrostatic pressure at the feed point where the local O 2 demand was highest.Unfortunately, the works referenced here did not include spatial profiles of DOT, and as such, direct comparison was not possible.O 2 limitations were not detected in any of the referenced works directly with the probes in the middle of the reactors.However, Bylund et al. (1998) hypothesized that O 2 limitations could have occurred around the feed points in their experiments, and likewise, Xu, Jahic, Blomsten et al. (1999) estimated based on the formate accumulation that approximately 12% of the culture volume would have been anoxic.The modeling performed here was in accordance with these hypotheses: a poorly mixing substrate feed might localize the O 2 demand such that the limitation is undetected by the electrode(s) just as suggested in the literature (Figure 5b).
1 were also used to estimate time-averaged efficiency factors for the experiments.Both the oxygen uptake and adaptation efficiency factors (Equations 18 and 20) were calculated for the same reported time points as in Figure 3c,d using the same substrate time-scale ratios calculated from the reported volumetric feed rates and substrate concentrations.In comparison with small scale, Xu, Jahic, Blomsten et al. (1999) reported a 0.31∕0.41≈ 76% yield efficiency on a large scale during the constant-feed phase.Time-averaged yield efficiencies of 88% and 86% were estimated here by oxygen uptake and adaptation efficiencies of 83% and 70%, respectively, resulting in a total yield effectivity of 0.88 × 0.86 = 75%.Bylund et al. (1998) Model estimations of axial distributions in the top-(t) and bottom-fed (b) large-scale cultivation experiments reported by Bylund ) and it also simplifies to standard Monod kinetics in homogeneous conditions.The linearization performed best in accordance with the Monod kinetics when the mean substrate concentration was at most of the same magnitude as the Monod constant (FiguresS1, S2, and S5).It has been pointed out previously that Monod kinetics have been validated in homogeneous conditions in chemostat and batch occur even with only t = 10s 95 mixing times characteristic of small-scale equipment.In large-scale reactors, where t > 200s 95 and longer mixing times are possible, mixing limitations may appear already with low biomass concentrations of 5gL −1 or feed rates of 1gL h −1 −1 (Tables Referenced large-scale experiments.
S−1 .The parameters and biomass concentrations necessary for calculating the rate constants were directly available in each study.T A B L E 1

10s 95 t = 100s 95 t = 200s
− ions at the controlled mean pH (Table1) were used as the initial state, and their local equilibrium concentrations were calculated such that Equation (27) applied to local CO 2 (aq) concentration as the equilibrium state.The mean pH and mean CO 2 (aq) concentra- tion were used to determine the mean HCO 3 − concentration.The local pH could then be obtained from either of the buffer equilibria.
Table S2 exemplifies the procedure.