Challenges in transfer of gas-liquid reactions from batch to continuous operation: dimensional analysis and simulations for aerobic oxidation

The transfer of gas-liquid reactions from conventional batch processes into continuous operation using milli and micro reactors is claimed as an important step towards process intensification. Importantly, this transfer step should be realized in an early phase of process development, already, in order to minimize research efforts towards the undesired operation strategy. The main challenge of this approach, therefore, arises from lack of knowledge in the early stage of process development and the resulting system with high degrees of freedom. This contribution presents an approach to tackle this challenge by means of mathematical modelling and simulation for the aerobic oxidation of 9,10-dihydroanthracene (DHA) catalyzed by polyoxometalates (POMs) being used as example for gas-liquid reactions. The reaction was chosen as it provides sufficient complexity, since it consists of three consecutive oxidation steps of DHA and a parallel catalytic redox-cycle according to a Mars-van-Krevelen mechanism. It also provides the challenge of unknown reaction kinetics, which have been estimated in this contribution. The dimensionless balance equations for reactor modeling are derived and parametrized based on early stage experimental results obtained in batch operation mode. The discrimination between batch and continuous operation was performed by means of characteristic dimensionless numbers using the identical mathematical model for comparability reasons. The model was used to perform sensitivity studies with emphasis on the interplay between mass transfer characteristics and reaction kinetics for both the batch and continuous operation mode. The simulation results show that the performance of both operation modes mainly depend on the oxidation state of the POM catalyst, which is caused by the differences in oxygen availability. Therefore, results obtained in batch operation mode are prone to be masked by mass transfer issues, which affects catalyst and reactor development at the same time and may thus cause maldevelopments. With respect to process development it can thus be concluded that the transfer from batch to continuous operation together with mathematical modeling is important in an early phase, already, in order to detect limitations misleading the development. Finally, even simple models with roughly estimated parameters from preliminary experiments are shown to be sufficient in the early phase and can systematically be improved, in the subsequent phases.


Introduction
Over the last two decades, milli and micro reactors receive an increasing attention due to the possibilities to shorten the way from lab to the market. Milli and micro reactors can be scaledup by numbering-up and the small characteristic lengths provide a good reaction control and thermal management [1]. For the production of fine chemicals and pharmaceutical products, the application of small-scale continuous reactors are already part of the reaction engineer's toolbox [2]. However, investigations in batch reactors are still the method of choice for highly viscous and sticky suspensions or for reaction mixtures probably producing precipitates, where easy cleaning is Article highlights • Batch-wise and continuously operated reactors are compared based on dimensionless reactor models and sensitivity analysis. • The dimensionless mass-transfer parameter is introduced to distinguish between reactors operated in batch and continuous mode. • Catalyst re-oxidation even limits slow Mars-van-Krevelen driven reactions, due to low oxygen availability.
favored over high reaction control [3]. Sometimes flow chemistry is claimed to be superior in general [4], although there are still challenges remaining, which can be circumvented by using batch reactors (e.g. chemical compatibility, clogging and fouling) [5]. Currently, liquid phase oxidations are under intensive investigation, because the standard processes for the conversion of petroleum-based raw materials will be substituted by bioderived raw materials in future [6]. Aerobic oxidation of organic substrates is particularly attractive, since it allows for functionalizing organic molecules and their implementation into a renewable chemical value chain. This typical application, however, comes along with hazardous risks if the reaction ignites, which becomes severe for large volume batch reactors [7]. One approach to minimize the risk of explosion or ignition is the utilization of the small characteristic diameters provided in milli or micro reactors, which can be chosen to be below the safe diameters of the explosive reaction mixture [8,9], however, microreactors are not inherently safe [9]. Therefore, hazardous operation windows are usually avoided. One of the most frequently used approaches is to operate below the limiting oxygen concentration (LOC) of the organic solvent and/or reagents used in a chemical reaction [10].
Despite these advantages for continuous operation in milli or micro reactors, the batch reactor is often chosen in the early stages of process development, due to the simpler handling and flexibility of operation conditions. Further reasons for preferring the batch against the flow reactor, especially for gas-liquid reactions, are the complexity and the realization of the desired flow pattern for realistic reaction mixtures. In addition, the specific mass transfer properties depend on many factors, such as physical properties of the fluids [11][12][13], the flow pattern itself and the reactor geometry, which induce significant uncertainties in the mass transfer correlations and predictability [14][15][16]. Nevertheless, flow reactors offer the opportunity to apply spectroscopic methods providing spatial concentration profiles, which is improving the knowledge gain during process development [17,18]. Furthermore, the understanding of fluid dynamics and mixing inside the milli and micro channels is improved, recently [19,20]. This progress together with the combination of computational fluid dynamics (CFD) and spectroscopic methods allows to analyze the mass transfer characteristics and thus minimize the uncertainties in prediction [21,22].
During development of gas-liquid reaction processes mathematical models play a key role for scale-up and reactor design (e.g. [23][24][25][26][27][28]), which need to be adapted to the specific case. Therefore, the design and operation parameters of the reactor, as well as the kinetic parameters of the reactions and mass transfer are required. Furthermore, a mass transfer model suitable for the problem has to be chosen [29]. In case of wellknown systems, the relevant information is available and the model can be validated by experimental data, before being used in reactor design and scale-up. Romanainen and Salmi worked on the simulation of tank reactors and bubble columns, where they implemented gas-liquid reactions applying film theory [23,24]. The focus was put on different strategies to solve the material balance numerically, which led to the conclusion, that a simultaneous calculation of bulk and film might be favorable, however, this strategic approach shows some instability for the numerical solution [23]. They also showed the application, demonstrating the simulation approach for the chlorination of butanoic acid in a bubble column. A penetration model was implemented by van Elk et al., who tried to overcome the asymptotic cases usually applied and demonstrated the possibility to cover more realistic application scenarios with a dynamic model [25]. Nevertheless, it is state of the art to apply simplified model equations or its analytical solutions to characterize a gas-liquid reaction system [30]. For instance, second order reactions [26] or the influence on reactant depletion in a batch reactor [27] are investigated based in the film model. However, Kucka et al. already derived expressions for the characteristics of the above mentioned models and the surface renewal model [28]. In their investigation, different methods for the determination of reaction kinetics are applied to experimental results from a stirred cell reactor and used to evaluate the reaction of CO 2 with Monoethanolamine (MEA). In the early stage of process development, however, the reaction kinetics might even be unknown and need to be estimated from preliminary experimental results. Markoš et al., for instance, present a very detailed dimensionless model based on the film theory for the oxidation of xylene. They point out that the agreement between simulation and experimental results depends on the quality of the used kinetics [31]. This is not surprising, since uncertainties are expected to affect simulation results from rigorous models. It has to be emphasized, however, that the benefit of modeling and simulation is not restricted to well-known systems, but also suitable for extrapolation and even for exploration towards unknown reaction systems. Modeling and simulation studies for those systems with high degree of freedom, though, is rare. Modelling and simulation for even more complex systems currently investigated, like three-phase systems, are far behind the experimental development [32,33]. A most recent investigation by Kappe and Coworkers is a positive example at the matured stage of investigation in the field of aerobic oxidation [34]. Nevertheless, the work shows the effort, which is necessary to keep astride experiments and simulation, resulting in a fit-for-purpose process model in an early stage of development. They combined reaction experiments with two-directional computational fluid dynamics (CFD) simulation for the selective oxidation of diphenyl sulfide to diphenyl sulfoxide applying a copper catalyst. The residence time distribution was characterized by a pulse experiment and the mass transfer coefficient was estimated using Higbie's penetration model. Thus, the investigation is combining the experimental investigations and simulation in order to scale and predict possible process conditions. Eventually, the coexistence of suitable models with experimental effort should be a characteristic commitment in reaction engineering.
The present contribution provides a sensitivity study based on modeling and simulation for the transfer of gas-liquid reactions from batch to continuous operation with emphasis on the particularities in an early phase of process development. Therefore, the dimensionless material balances are derived and parametrized based on experimental results obtained in batch operation mode. The discrimination between batch and continuous operation was performed by means of characteristic dimensionless numbers using the identical mathematical model for comparability reasons. The aerobic oxidation of 9,10-dihydroanthracene (DHA) catalyzed by polyoxometalates (POMs) is used as an example reaction, being subject of recent research efforts in catalyst [35,36] and reactor [37] development. The chosen example provides sufficient complexity in order to demonstrate the potential of the model-driven batch-to-conti approach in early stages of process development, since an organic reactant undergoes consecutive oxidation steps catalyzed by POMs via a Marsvan-Krevelen (MvK) mechanism. In addition, neither the reaction mechanism, nor the kinetics are known in sufficient detail, which is usually the case, if catalyst development is in the focus of the research.
We show that a model-based exploration of the complex gas-liquid system is possible, despite the presence of uncertainties. In particular, we observed that the oxidation state of the POM catalyst is the crucial factor in the transfer from batch to continuous operation, which is caused by the differences in oxygen availability. Importantly, this observation is not easily accessible by experiments, only, but requires support by modeling and simulation. We therefore demonstrate that even simple models with roughly estimated parameters from preliminary experiments are very important tools for a simulation-driven transfer from batch to continuous operation already in an early stage of process development.

General concept
In order to describe the "batch-to-conti" transformation on a comparable basis the following approach is used to set up the balance equations. In particular, both the stirred tank reactor operated in batch mode (bSTR) and the continuously operated plug flow reactor (cPFR) are modeled based on the equal meaning of reaction time t in the batch reactor and residence time τ in the plug flow reactor. Furthermore, the reactor model consists of the combination of balance equations at two length scales: The meso scale refers to the balances in proximity of the phase boundary, while the macro or reactor scale covers the balances for the fluid bulk.
The example reaction chosen justifies the assumption of isothermal conditions, since the adiabatic temperature rise of the reaction mixture is small (21 K) and therefore, heat balances are neglected. Furthermore, film theory is assumed to describe the transport processes in proximity of the phase boundary, in order to allow for illustrative presentation of the simulation results. Concerning the temporal behavior, steady-state conditions are assumed for the meso scale balances. In contrast, unsteady-state balances are formulated for the macro scale, which allows to utilize the residence time to be used as comparison basis for both reactor types.

Balance equations
Based on the main assumptions, the material balance at the meso scale can be formulated according to Eq. (1) for a flat geometry of the liquid film with a thickness of δ L , in order to derive the profile of the concentration of component A i in the liquid film c i,L as function of the spatial coordinate x. The rate r j,L of reaction j in the liquid film is a function of the spatial coordinate x, as well, as it depends on the concentrations of the reactants. Furthermore, D i,L refers to the diffusion coefficient in the liquid phase, while ν i,j holds for the stoichiometric coefficients.
The boundary condition at the gas-liquid interface (x = 0) is provided in Eq. (3), which assumes no chemical reaction taking place in the gas film. Here, k i,G refers to the gas sided mass transfer coefficient and c * i;L to the concentration at the gasliquid interface. Furthermore, Henry's law is used to convert the partial pressure of component A i in the gas phase p i,G into the respective equilibrium concentration c i,G,eq with the Henry-coefficient H i according to Eq. (2). For components with negligible vapor pressures (H i = 0) the boundary condition can be simplified to Eq. (4) leading to the absence of mass transfer across the gas-liquid interface. At the interface between the liquid film and the liquid bulk (x = δ L ) the concentration has to be a continuous function and thus Eq. (5) is used as respective boundary condition. Note that the bulk fluid is assumed to be well mixed in direction of the coordinate x according to the film theory.
At the macro scale Eqs. (6) and (7) describe the material balances of the gas and liquid bulk, respectively. Note that the left side of both equations refer to the accumulation in the bulk fluid. In case of the bSTR τ refers to the time dimension and thus the reaction time, while it corresponds to the residence time τ = z/u TP with z being the spatial coordinate in flow direction and u TP the two-phase velocity of gas and liquid in the cPFR. Therefore, the variable τ quantifies the reaction time for both the bSTR and cPFR, which allows for direct comparison of the results obtained for both reactors with respect to reaction time.
The right side of Eq. (6) describes the mass transfer between the gas bulk and the gas-liquid interface according to the boundary condition introduced in Eq. (3) with the specific gas-liquid interfacial surface area a G,L = A G,L /V R , the gas holdup ε G = V G /V R and the reactor volume V R . In Eq. (7) chemical reaction in the liquid bulk r j,L,b is considered, which equals r j,L (x = δ L ). Furthermore, J i,δ refers to the mass flux density of component A i between the liquid film and the liquid bulk according to Eq. (8). The balance equations at the macro scale thus correspond to an initial value problem with the initial conditions provided in Eqs. (9) and (10).

Dimensionless balance equations
The dimensionless form of the balance equations is derived based on the following definitions. First, the concentrations and initial conditions are transformed into residual fractions f i according to Eqs. (11) and (12) with c 1,ref chosen as the reference concentration. Furthermore, the dimensionless spatial coordinate χ in the film is defined as indicated in Eq. (13).
The dimensionless time θ is obtained as the ratio between the reaction time τ and the overall reaction time τ in both reactor types (Eq. (14)). The dimensionless reaction rate ω j is defined with respect to a reference reaction rate r ref according to Eq. (15). It has to be mentioned that the reference reaction rate can arbitrarily be chosen and corresponds to the kinetics of the virtual reaction as stated in Eq. (31) for the present contribution, as will be discussed below in detail.
Furthermore, the Hatta number Ha and the Damköhler number Da I are defined (Eqs. (16) and (17)). It becomes obvious that the chosen reference reaction rate determines the meaning of the values obtained for Ha and Da I and thus requires a closer look as discussed in detail below (see Eq. (31)). For Ha the assumption of film theory with k i,L = D i,L /δ L allows to eliminate the liquid film thickness δ L . The Biot number Bi i relates the gas and liquid sided mass transfer to each other (Eq. (18)). Finally, the Hinterland ratio Hi relates the liquid volume V L to the volume of the liquid film V L,f and can be expressed using the film theory and geometric considerations as shown in Eq. (19). Hi The resulting dimensionless material balance at the meso scale is shown in Eq. (20) with the corresponding boundary conditions in Eqs. (21) and (22). Note that steady-state is assumed for these balance equations.
At the macro scale Eqs. (23) and (24) are derived accordingly with the initial conditions provided in Eqs. (25) and (26). For further simplification, we assume that the gas phase consists of one single component only, which leads to constant value for f i,G as indicated in Eq. (23). Finally, we assume an ideal segmented flow pattern for the gas-liquid two-phase flow in the cPFR (see Fig. 1) with well-mixed gas and liquid slugs and neglecting the film between the gas slug and the wall. Note that one unit cell can be regarded as control volume for the segmented flow model assuming the unit cell length being much smaller than the reactor length.
Implementation and simulation

Challenges from reaction network and catalysis
The example reaction chosen for this study is the oxidation of 9,10-dihydroanthracene (C 14 H 12 ) to anthraquinone (C 14 H 8 O 2 ) via the intermediates anthracene (C 14 H 10 ) and anthrone (C 14 H 10 O) using gaseous oxygen as reactant according to the stoichiometric Eqs. (R1) to (R3) [36]. The by-product water is sufficiently soluble in N,Ndimethylformamid (DMF) used as solvent, in which the chemical reaction takes place. All reaction steps are assumed to be catalyzed by a POM ({Ba 4 V 14 O 38 }), which provides oxygen to the organic substrate and is re-oxidized by molecular oxygen based on a MvK mechanism (Eq. (R4)).
This reaction network is highly relevant, due to the importance of anthraquinone in the chemical value chain, being currently produced by cycloaddition or coupling reactions, e.g. by Friedel-Crafts-reaction of benzene and phthalic anhydride. From theoretical point of view and more relevant to the current study, the reaction network provides several interesting features. Importantly, the consecutive character is well understood in chemical reaction engineering, which provides the basis in order to perform fundamental studies. One aspect arises from the POM being present in the liquid phase in form of clusters in the size of 3025 g mol −1 , which can therefore be considered as a molecular, rather than a solid catalyst. This renders the system being homogeneously catalyzed, for which the theoretical framework of non-catalytic gas-liquid reactions can be applied counterintuitively. In particular, the implementation of the MvK mechanism in modeling of gas-liquid reactions is not yet performed in scientific literature to the best of our knowledge.
Another interesting aspect is motivated by the typical challenge in research for derivation of kinetic data, where experiments have to be designed despite a significant lack of knowledge on the reaction system, especially in an early stage of research. For multi-phase reactions, in particular, care has to be taken in order to avoid the obtained kinetic data to be falsified by heat and mass transfer issues. Even though several rules of thumb are available to tackle this challenge (e.g. the Weisz-Prater criterion in heterogeneous catalysis), these rules are usually applicable to very simplified reaction systems only. In this context the reaction system chosen is too complex and still not fully understood, neither mechanistically, nor with respect to the reaction kinetics. On top of this, the variation of the POM based catalyst during research progress may affect both the reaction mechanism and the kinetics, which hampers the transfer of knowledge gained so far.
The scientific key question for experimental derivation of reliable kinetic data are the concentration and temperature Fig. 1 Model for the unit cell in a segmented flow (Taylor-flow) regime, with the unit cell consisting of one gas bubble and one liquid slug, as well as an illustration of the two film model at the meso scale profiles in proximity of the gas-liquid interface, since they allow to distinguish reaction kinetics from the effects of transport limitations. Operando techniques are not capable to solve these issues completely yet [22,38], since resolution of those profiles is still a matter of research [39,40]. Therefore, modeling and simulation is a promising approach, in order to support the design of experiments and interpretation of the results, as will be demonstrated in the present contribution.

Reaction system
For implementation, the matrix of stoichiometric coefficients is derived from the stoichiometric reaction Eqs. (R1) to (R4) as summarized in Table 1. In addition, the components A i are also numbered in the following for simplification of the notation. Note that oxygen is present in the gas phase and also soluble in the liquid, while all other components are considered to be present in the liquid phase only, due to a negligible vapor pressure of the organic reactants and the catalyst.
The kinetics involved for the occurring reactions are not known yet and thus require to be estimated. Therefore, we chose a simple approach assuming all reactions to be irreversible and first order in the respective reactants (Table 2). In particular, the catalyst is considered to take part in each reaction, either as oxygen donor or acceptor according to the MvK mechanism.
Reaction (R4) is chosen as the reference, being first order both in oxygen and catalyst concentration according to Eq. (27). Therefore, the initial equilibrium concentration of molecular oxygen in the gas bulk corresponds to the reference concentration c 1,ref , which is also used for definition of the dimensionless numbers and variables above. Furthermore, the total catalyst concentration, comprising the oxidized and reduced state, is considered as c cat,ref , which is constant throughout the reaction progress. The kinetic constant of the reference reaction k ref needs to be scaled for sensitivity analysis during simulation experiments or to meet experimental data. Note that the effect of reaction temperature is included in k ref , only, for the presented isothermal balance equations.
We have chosen this type of reference reaction kinetics, as it involves the gaseous reactant on the one hand, which has to be transported across the gas-liquid interface. On the other hand, it also includes the concentration of a reactant fed with the liquid phase, being the catalyst in this case. The reference kinetics thus put the most important reactants of the complex reaction network in the center of the considerations, being the catalyst in the focus of ongoing research. It has to be mentioned again that the choice of the reference reaction kinetics depends on the individual case and can even describe a virtual reaction not taking place in real experiments, since it is only required to derive the dimensionless balance equations. After all, the meaning of the specific values for the dimensionless numbers is defined by the choice of the reference.

Parametrization
For parametrization the dimensionless numbers Da I and Ha are estimated based on experiments performed in the bSTR by Lechner et al. [36]. The overall reaction time in those experiments was τ ¼ 7 h, with an initial oxygen partial pressure of p O 2 ;0 ¼ 8 bar as well as the initial concentration of component A 2 of c 2,L,b,0 = 45.8 mol m −3 and of the catalyst of c 7,L,b,0 = 0.9 mol m −3 in the liquid phase. Components A 3 to A 6 and A 8 Table 1 Matrix of stoichiometric coefficients of the reaction system (the solvent DMF is implemented as inert component A 9 ) are not present in the reactor initially, while the total liquid volume was V L = 5 · 10 −5 m 3 . The authors report an initial conversion of component A 2 obtained after a reaction time of Δτ = 1 h of Δc 2,L,b = 11.45 mol m −3 , as well as a formation of component A 4 of Δc 4,L,b = 5.95 mol m −3 and A 5 of Δc 5,L,b = 5.50 mol m −3 . From these changes in concentration the number of catalytic cycles of the catalyst, which is the turn-over number (TON), can be calculated according to Eq. (28). Note, that the first term in the numerator represents all turn-overs in the first reaction.
With the Henry coefficient of [41,42], as well as the estimated orders of magnitudes for the diffusion coefficient of O 2 in the liquid phase of D O 2 ;L ≈ 5 Á 10 −9 m 2 s −1 [43] and for the liquid film thickness of δ L ≈ 40 μm, Table 3 summarizes the relevant values for the dimensionless numbers of a typical experiment in a bSTR. The diffusivity of all other components is assumed to be five times smaller, compared to oxygen diffusivity, however, diffusivities for the catalyst are unavailable and might be even smaller than for the Anthracene derivates (D A 3 ;L ¼ D A 5 ;L ≈ 1 Á 10 −9 m 2 s −1 [44,45]). Furthermore, the gas-sided mass transfer is assumed to be not limiting (Bi i → ∞).

Reactor comparison
Interestingly, the dimensionless balance equations allow to express the difference between bSTR and cPFR by the value for Hi only, since Eq. (19) can be rewritten to Eq. (29) assuming film theory. It becomes obvious that Hi consists of properties strongly depending on the reactor type used and is thus suitable to distinguish between batch or continuous operation. In particular, Hi becomes small for segmented flow in the cPFR, due to the high values for specific interfacial surface area. The estimation of the order of magnitude for ε L ≈ 0.5 and a G,L ≈ 1000 m −1 typical for segmented flow (including liquid film surrounding the bubbles) in channels of 1.6 mm (1/16 in.) diameter, as well as for δ L ≈ 40 μm (rough estimate from k L a G,L values from [15]), gives Hi ≈ 12.5. In contrast, a Hinterland ratio of Hi ≈ 250 is obtained for the used bSTR with ε L ≈ 0.5, a G,L ≈ 50 m −1 and δ L ≈ 40 μm (agitated tank with small energy input, c.f. [46]).
Note that Da I is not suitable to distinguish between the operation strategies, as it relates the reaction time with the intrinsic reaction kinetics, which is identical for both reactors. In contrast, Ha quantifies the mass transfer characteristics and may thus weakly depend on the reactor type chosen, only.
The interplay between reaction in the bulk and the profiles in the film can illustratively be discussed based on the mass transport parameter T (Eq. (30)). Comparison with Eq. (24) exhibits that this parameter describes the transport of the limiting component supplied by the gas phase A 1 from the liquid film into the liquid bulk and therefore links the meso and macroscale with each other. It also becomes obvious from Eq. (24) that for high values of T the gradient of the residual fraction at χ = 1 can be small for constant flux densities, which results in flat profiles for A 1 in the film.
This parameter is also suitable for fundamental studies on the batch-to-conti approach, since it allows a direct comparison of results obtained for various values of Hi, Ha, and Da I on a common basis by lumping the individual effects. In other words, since all these parameters are affected by changing the operation strategy from batch to continuous operation, the proper comparison would require to perform an intensive study with these three parameters. Lumping the individual effects in T reduces the parameter space without losing significance.

Implementation and simulation
The model equations are implemented to Matlab (version R2019b) and solved with ode15s (4th order backward differentiation formula -BDF) in order to perform a sensitivity analysis by simulation experiments. The simulations focus on the comparison between bSTR and cPFR, which are distinguished based on Hi, while keeping all other parameters equal for comparability reasons. Furthermore, the initial conditions used in the above-mentioned experiment are applied for simulations, as well. The diffusion coefficients of all Anthracene derivates and the catalyst in the liquid phase are assumed to be equal for simplicity reasons, though this assumption may induce errors in the obtained results. We therefore focus the discussion of the results on the trends observed, rather than on the interpretation of specific numbers. For evaluation of the reactor performance, following definitions are used for conversion X of component A 2 and the yield Y i towards component A i (Eq. (31)). The degree of reduction of the catalyst α is quantified by Eq. (32).
Results and discussion

Base case
The base case assumes that the dimensionless initial rates of the reactions (R1) to (R3) amount to ω j,0 = 0.1 for j ∈ {1, 2, 3}, which means that those reactions are slower than the reference reaction (R4). This scenario corresponds to the typical approximation used in homogeneous catalysis, which neglects the formation of concentration gradients of the molecular catalyst. Indeed, it is known from experiments, that the specific chemistry of catalytic steps is a delicate task (c.f. [19]), but is typically neglected. In the present example, however, the catalyst is stoichiometrically involved in the reaction, due to the MvK mechanism proposed. Therefore, certainly pronounced gradients in oxidized and reduced catalyst species may evolve, depending on the reaction rates.
The assumption limits those gradients on the one hand, but offers the respective discussion at the same time. Therefore, this approach allows to shed light on the importance of the spatial distribution of the catalyst oxidation state within the reactor on the measured results. In order to fulfill the assumption, the kinetic c o n s t a n t s r e q u i r e t o a m o u n t t o k j = 2 . 2 8 · 10 −5 m 3 mol −1 s −1 for j ∈ {1, 2, 3}. Figure 2 shows the simulation results for the bSTR under base case conditions. The observed profiles of the residual fractions for the organic reactants as function of reaction time (Fig. 2a) are typical for the consecutive reaction considered in this example and agree qualitatively with experimental results [36]. After 7 h of reaction (θ = 1) the conversion of A 2 reaches ca. 88%, while the yield for all anthracene derivates is comparable in the range of 25 to 35%. The high numbers for the residual fractions are caused by the low solubility of O 2 in DMF and thus a small value for c 1,ref . One may conclude that the solubility of O 2 appears to be one important factor responsible for the long reaction times required in the experiments in order to achieve reasonable conversions. In Fig. 2c the corresponding profiles of the oxidized and reduced catalyst species are shown. The initially fully oxidized catalyst is reduced rapidly in the very early stage of the reaction and remains at a rather constant residual fraction of the oxidized state of f 7,L,b ≈ 16 thereafter. The profiles of the residual fraction of oxygen and the degree of reduction of the catalyst α at the meso scale are shown in Fig. 2b for three different reaction times θ I , θ II and θ III indicated in Fig. 2a (corresponding to ca. 30 min, 3.5 h and 7 h). Obviously, the degree of reduction of the catalyst is constant in the liquid film, as well. The oxygen profile, however, declines almost linearly from the phase boundary at χ = 0 towards the liquid bulk at χ = 1, which agrees with the expectations from the small Ha used for the simulations.
The almost constant degree of reduction of the catalyst α at meso and macro scale renders the apparent kinetics for (R1) to (R3) being first order in concentration of the respective organic reactant. The oxygen availability in the liquid phase, however, is the limiting factor for the re-oxidation of the catalyst in (R4) and hence for the overall reaction process. In particular, the residual fraction of oxygen in the liquid bulk amounts to ca. 0.2 only. Considering that the reaction mainly takes place in the bulk for the large Hi used in the simulations, the limiting effect of oxygen availability is not surprising. Importantly, the catalyst is provided in large excess with respect to oxygen (κ 7,L = 22.9) and therefore the degree of reduction α is rather independent on the reaction progress and limited to ca. 0.3. The negligible differences in the profiles observed at the meso scale are caused by the slight changes of α with reaction time at the macro scale.
The simulation results for the base case in the cPFR are shown in Fig. 3, providing qualitatively similar profiles compared to the bSTR with a slightly higher conversion for A 2 of ca. 95% and yields in the range of 14 to 59% for all organic products. The degree of reduction of the catalyst is also almost constant in the film (Fig. 3b) and with reaction time (Fig. 3c), but lower than in the bSTR. Furthermore, the profiles at the meso scale are not affected by the reaction time. Compared to the bSTR, however, the profile of oxygen residual fraction in the film is almost constant and close to unity (ca. 95%) even in the liquid bulk. Hence, the oxygen availability is hardly limiting, which consequently means that the re-oxidation rate of the catalyst in the bulk can be regarded as the limiting factor. The differently pronounced oxygen profiles in the film in the base case simulations can be explained by the parameter T, which is significantly smaller for the bSTR (T bSTR = 320) compared to the CFPR (T cPFR = 6400). It has to be mentioned that the convection and reaction terms in Eq. (24) are rather unaffected by the reactor choice given by the similar profiles of the residual fractions with reaction time. Therefore, the differences in the oxygen profile in the film can directly be attributed to the mass transport parameter T.
The degree of catalyst reduction α can be predicted depending on the reaction extent assuming that the rate of catalyst reduction in reactions (R1) to (R3) equals the rate of reoxidation in (R4). This assumption, expressed in Eq. (35), requires negligible changes in degree of reduction with spatial coordinate and reaction time. After introduction of the reaction kinetics (Table 2) and the definition of the reduction degree (Eq. (32)) Eqs. (36) to (38) can be derived, respectively.
Equation (38) provides the basis for an estimate of the expected degree of reduction, in order to design and interpret experimental results. Most importantly, the relevance of oxygen solubility and re-oxidation rate becomes apparent both expressed by the term k ref c 1,G,eq , since increasing this term leads to a decrease in the degree of reduction, which is beneficial for the reactor performance. Even though this conclusion is drawn based on the specific reaction kinetics assumed in the present example, it is of general importance for similar reaction systems and other kinetics, as well. Eventually, the relationship shows that in-situ measurements of the catalyst state will give a further insight into the underlying kinetics.

Kinetics
In order to represent the experimental results reported by Lechner et al. [36] the dimensionless reaction rates have been adjusted according to the following values: ω 1,0 = 0.5, ω 2,0 = 25, and ω 3,0 = 5. Note that these values are not obtained by fitting to experimental data, due to lack of knowledge in reaction kinetics. Instead, emphasis is on studying the order of magnitude effects of the reaction rates on the observations. The chosen values mean that the rate of reaction (R1) is limiting, while the selectivity towards the intermediate (A 3 and A 4 ) and final (A 5 ) products are governed by the respective dimensionless reaction rates. In particular, A 3 formation is expected to be hardly detectable, since it is directly converted into A 4 given by the high ω 2,0 /ω 1,0 ratio.
The direct comparison of experimental and simulation results obtained for both the bSTR and cPFR reactor is summarized in Table 4 and depicted in Fig. 4. Note that the simulations are based on identical parametrization, except for Hi for both reactor types. Therefore, the direct comparison based on simulations is possible, since the errors induced by unknown kinetics and further simplifying assumptions contribute equally to the obtained results for both reactors. From Table 4 it can be observed that the simulation is capable of predicting the trend of the experimental data, while the specific values deviate to a certain extent. The discrepancy can be explained by the fact that the reaction kinetics are not known and therefore approximated by simple power-law approaches. On the other hand, this simple approach obviously allows to gain meaningful insights into the reaction system and the reactor behavior. Figure 4 provides the direct comparison of both reactor types based on the profiles of the organic reactants with reaction time obtained by simulations. The observed trends are comparable for bSTR and cPFR, exhibiting a declining residual fraction for A 2 , negligible values for A 3 , and monotonously increasing values for A 5 . Component A 4 exhibits a profile typical for intermediates in consecutive reactions. The resulting profiles in residual fractions thus agree with the expectations associated with the chosen values for the dimensionless reaction rates in a consecutive reaction network.
The detailed comparison reveals that the cPFR exhibits better performance with respect to conversion of A 2 and yield of A 5 obtained at reaction time of θ = 1. Comparing the reaction time required to reach a conversion of X = 0.2 and the yield in A 5 at that reaction time, the cPFR (θ cPFR = 0.15, Y 5,cPFR = 0.11) outperforms the bSTR (θ bSTR = 0.25, Table 4 Comparison of reactor performance in simulation and experiment [36] for both reactor types (no cPFR experiments with fluidic residence time τ fl > 1 h available)

Simulation Experiment
Reaction time Reactor  Y 5,bSTR = 0.11) clearly. The reason for this observation has to be related to the value of Hi, since this parameter distinguishes between both reactor types, only. In particular, Hi is responsible for oxygen availability in the liquid bulk and hence for the degree of reduction of the catalyst, as discussed for the base case. Therefore, the profiles of the oxygen residual fraction together with the degree of reduction as function of reaction time are depicted in Fig. 5 for cPFR and bSTR. Obviously, the degree of reduction is significantly lower in the cPFR reaching a maximum value of about 15%, while that in the bSTR even exceeds 50%. This observation is associated with the oxygen residual fraction, which is rather constant with reaction time at ca. 0.9 and 0.2 for the cPFR and bSTR, respectively. Interestingly, the oxygen residual fraction is rather constant, while the degree of reduction exhibit significant changes with reaction time. This can be explained by the high value for ω 2,0 , which means that the catalyst reduction in reaction (R2) is twenty-five times faster than the catalyst reoxidation in reaction (R4) under initial conditions. Therefore, the oxygen consuming re-oxidation reaction is limiting in the simulated cases and the profile of reduction degree is governed by the rates of reactions (R1) to (R3). On the other hand, the oxygen transport through the film is limiting due to the high value for Hi in the bSTR, which leads to pronounced profiles in oxygen residual fraction in the film and thus low values in the bulk. This limits the re-oxidation, as well, which is first order in oxygen. Finally, it has to be considered that the oxygen demand depends on the reaction extent, due to differences in stoichiometry and the intrinsic rate of reactions (R1) to (R3) and thus changes with reaction time (e.g. Eq. (38)). The interplay of these factors is therefore responsible for the observed profiles and renders the cPFR beneficial for kinetic experiments, because mass transport limitations are less pronounced. Interestingly, the simulation results exhibiting a higher oxygen saturation and smaller degree of reduction in the cPFR agree in principle with experimental observations. The POM catalyst used in the experiments changes color from orange to green upon increasing degree of reduction and even becomes colorless after decomposition, which justifies the predictions from simulations qualitatively.

Sensitivity analysis in mass transfer
The effect of mass transfer is studied based on the lumped mass transport parameter T introduced in Eq. (30). For sensitivity analysis this parameter is varied between 30 and 120′ 000, while each value is justified by typical conditions for either bSTR or cPFR. As basis for a meaningful range of T the film thickness was varied between 2 and 400 μm, which even covers values beyond the expected boundaries for bSTR and typical flow regimes in cPFR. We thereby account for uncertainties in estimation of T by broadening of the considered value range and avoid the complexity arising from more accurate determination at the same time. Figure 6 depicts the impact of the mass transport parameter T on the maximum degree of catalyst reduction in the bulk achieved within the reaction time and the conversion of component A 2 after 55 min (θ = 0.13) and 7 h (θ = 1) of reaction. The solid lines correspond to typical conditions for each reactor type, while the dashed lines indicate the extrapolation towards unusual conditions, in order to support the direct comparison of batch and continuous operation. For illustration, the range of T corresponding to the film thickness between 2 and 400 μm is indicated for the bSTR and the cPFR, as well. It becomes apparent that both reactor types overlap in the range of 250 < T < 2000, which correspond to extrapolated values for each and therefore unlikely operation conditions. The degree of maximum reduction appears to be a function of the mass transport parameter T only, but does not depend on the reactor type, since the α max (T) profiles of cPFR and bSTR are clearly overlapping. Furthermore, it appears that extrapolation over the whole range of T leads to identical profiles irrespective of the reactor type. It can also be observed that the reduction degree is constant for the typical operation range in the cPFR, but exhibits a strong sensitivity towards the mass transport parameter for the bSTR. Importantly, the maximum reduction degree is very small and the average value therefore even smaller for the cPFR. Consequently, the catalyst can be assumed to be homogeneously distributed in the liquid phase in the oxidized state, which corresponds to the typical assumption for molecular catalysts in homogeneous catalysis. In contrast, the degree of reduction reaches very high values for the bSTR, which lead to the irreversible deactivation of the POM-based catalyst by decomposition. Therefore, gradients in the concentration of oxidized and reduced catalyst species probably develop. The reason for the high degrees of reduction is the limited availability of oxygen for re-oxidation for low values of T and thus reasonable mass transfer limitations for oxygen. Hence, the cPFR is beneficial for kinetic measurements, due to absence of severe mass transfer limitations at high T values. Considering the conversion, the results for the bSTR raise with increasing T and asymptotically reach those of the cPFR and therefore intrinsic conditions. Furthermore, the sensitivity towards the mass transport parameter is more expressed in the bSTR.
The sensitivity of both the conversion and degree of reduction in the bSTR indicate that a correlation between both exists: with decreasing degree of reduction the conversion increases. This is most probably caused by the enhanced oxygen availability for higher values of the mass transport parameter and thus improved re-oxidation of the catalyst by reaction (R4), which is first order in O 2 . This leads to a higher concentration of the oxidized catalyst species A 7 and thus an increased rate for the oxidation of A 2 via reaction (R1), which is first order in A 2 and A 7 . The rates of the consecutive reactions (R2) and (R3) profit from the increasing concentration of A 7 , as well, but form the reduced catalyst species stoichiometrically at the same time. For additional information, the reader is referred to the supporting information, where the degree of reduction and the rates are plotted for the highest and lowest value of the mass transport parameter, respectively.
Interestingly, the conversion in the cPFR is constant in the typical operation window, while it increases for smaller T values. The dimensionless formulation of the balance equations hampers a clear argumentation of this unexpected observation. Importantly, the mass transport parameter T results from lumping of the degrees of freedom of the reactor operation expressed by Hi, Ha and Da I , which allows different sets of those degrees of freedom to yield identical T values. For instance, a reduction of Hi from 250 to 12.5 for bSTR to cPFR can be compensated by a decrease in Da I from 10 to 0.5 in order to maintain T/Ha 2 = 400. The relation T/Da I , representing mass transfer into and the reaction within liquid bulk, changes from T/Da I = 0.4 in the first to T/Da I = 8 in the second case. This means that the interplay between reaction in the bulk and the processes in the liquid film play an important role for the performance of the chosen reactor. Therefore, fundamental differences between both cases exist depending on the specific parametrization, though the balance equations are identical. Along this argumentation, an increase in Ha leads to a decrease in T at constant Da I and Hi, which is associated with an increasing enhancement factor E and therefore an improved mass flux of oxygen across the gas-liquid interface. In the same vein, a reduction in T induces the gradient of the residual fractions at χ = 1 to raise (see Eq. (24)), which enhances the mass flux into the bulk, as well. According to Eq. (30) small T values can be achieved for fast reactions, expressed by high values for Ha and small values for Da I , and cause gradients at the meso scale to develop. Therefore, large values of T allow to neglect gradients in the film and thus to assume intrinsic conditions, while decreasing T values lead to deviation from intrinsic operating points. The latter situation thus requires more sophisticated models for analysis of experimental data, while a one-dimensional, pseudo-homogeneous model is sufficient for the intrinsic case. It also becomes apparent that a straight-forward argumentation of the obtained results based on the dimensionless numbers only, is misleading, due to the interdependencies of their definitions.

Variation of kinetics in the reaction network
The impact of the kinetics in the consecutive reaction network is studied in two scenarios. In scenario I the kinetics of reaction (R1) is varied according to ω 1,0 = ∈ [0.005, 5], while the kinetics of reaction (R2) and (R3) are kept constant at ω 2,0 = 25 and ω 3,0 = 5. This scenario therefore covers various cases relating the reduction of the catalyst in reaction (R1) to its reoxidation in reaction (R4). Within the consecutive network Fig. 6 Maximum degree of catalyst reduction α max and conversion in A 2 in the liquid bulk as function of the mass transport parameter T for the bSTR (squares) and cPFR (circles) with typical (solid lines, filled symbols) and extrapolated (dashed lines, open symbols) operation range; lines are shown to guide the eye, typical range in film thickness δ L related to mass transport parameter T is indicated by bars for both reactor types reaction (R1) is thus always limiting, since reaction (R2) is at least five times faster. It has to be considered, however, that reaction (R2) and (R3) require oxygen stoichiometrically, as well, and thus contribute to catalyst reduction.
The obtained simulation results (Fig. 7a) exhibit similar and almost negligible degrees of reduction for bSTR and cPFR for very small values of ω 1,0 associated with similar values for conversion and yield. Therefore, oxygen availability appears to be not limiting, neither in the molecular form required for the catalyst re-oxidation, nor in form of the oxidized catalyst needed for reaction (R1) to (R3), which renders the reaction rates being intrinsic. With increasing values for ω 1,0 the degree of reduction raises for both the bSTR and the cPFR, while the increase is significantly more pronounced for the bSTR. This can be explained by the increasing demand of oxygen in reaction (R1) to (R3) induced by the higher reaction rates, leading to the re-oxidation of the catalyst becoming limiting. The difference between both reactors is caused by the mass transfer parameter T, which is significantly smaller for the bSTR as discussed above. This leads to more pronounced mass transfer limitations and therefore to a lower oxygen availability. In addition to the observations at the macro scale shown in Fig. 7a gradients develop in the liquid film as discussed above for the base case.
The conversion X and yield Y 5 raises with ω 1,0 for the cPFR and reaches unity asymptotically. Furthermore, both values are very similar indicating a high selectivity towards component A 5 after the reaction time of θ = 1. In contrast, the conversion raises, while the yield exhibits a maximum for the bSTR, which means that the selectivity towards the intermediate components A 3 and A 4 increases with ω 1,0 . The observed maximum and subsequent decline in yield is associated with an increasing trend of the degree of catalyst reduction beyond α max > 0.6. Hence, in comparison to the cPFR the yield in the bSTR is obviously limited by the re-oxidation of the catalyst. In particular, the available oxygen is mainly consumed by reaction (R1) and (R2), since ω 2,0 /ω 3,0 = 5, and thus reaction (R3) suffers from oxygen depletion, which becomes more severe for high ω 1,0 values. Scenario II extends the previous analysis by variation of the rates of reaction (R1) to (R3) proportionally, according to ω 1,0 = 0.5 m, ω 2,0 = 25 m and ω 3,0 = 5 m with the activity ratio m = ∈ [0. 1,10]. It therefore scales all reactions contributing to catalyst reduction equally and thus relates the overall rates of catalyst reduction and re-oxidation to each other. Hence, the factor m expresses the activity ratio of the reduction and the reoxidation half-cycle during the MvK mechanism. The direct comparison to scenario I is possible for m = 1 and indicated in Fig. 7 as dashed vertical line.
The results in Fig. 7b exhibit an increasing degree of reduction with raising activity ratio, since the re-oxidation halfcycle of the catalyst becomes limiting. Furthermore, both the bSTR and the cPFR exhibit similar values of catalyst reduction degree, conversion and yield for m < 0.25, while the deviation becomes more pronounced for raising m. This can again by explained by the smaller mass transport parameter T in the bSTR, which limits the oxygen availability required for catalyst re-oxidation. Importantly, higher m values mean that catalyst re-oxidation is hampered under reference conditions at maximum oxygen availability already. During reaction oxygen is consumed and thus the catalyst reduction degree becomes even more sensitive towards mass transport restrictions. Interestingly, the selectivity towards component A 5 is high for both reactors expressed by the slight deviation between conversion and yield, since all reactions suffer proportionally from the higher reduction degree of the catalyst. The lower selectivity as well as the deviation of the maximum value for conversion and yield observed in the bSTR from the expected value of one for irreversible reactions indicates that the catalyst re-oxidation limits the reactor performance in that case.
The results of the kinetic study exhibit a very similar behavior of the bSTR and the cPFR in the intrinsic regime, as expected, which is beneficial for experimental investigation of kinetic parameters or catalyst development. The low values for conversion and yield can easily be compensated by increasing residence time in this regime. Under non-intrinsic conditions, however, pronounced gradients and oxygen transport limitations develop, which lead to several challenges in evaluation of experimental data. In case that catalyst reoxidation is limiting, for instance, the experimental determination of conversion and yield of the organic reactants becomes rather meaningless. Instead, emphasis should be put on the additional determination of the catalyst reduction degree, in order to evaluate the information obtained from the product composition data (see Eq. (38)).
Catalyst development often aims at increasing the activity with respect to conversion of the organic reactants, which is expressed by increasing ω 1,0 in scenario I. Interestingly, the results clearly show that an increase in activity is not necessarily correlated with the measured conversion or yield, which in particular holds for the non-intrinsic case with severe mass transport limitations. Obviously, the choice of an appropriate reactor and operations conditions become very important to provide a correlation between experimental results and development aims. The results in Fig. 7a, obtained for catalyst activity ranging over three orders of magnitude, clearly underline the benefits of the cPFR for that purpose.
From Fig. 7b it can be concluded that catalyst development requires a holistic approach, emphasizing both the reduction and re-oxidation half cycle. In other words, improving the activity towards conversion of the organic reactants require a proportional increase in catalyst re-oxidation activity, as well. Otherwise, the re-oxidation step becomes limiting and thus the catalyst reduction degree increases towards rather high values in the non-intrinsic case. This leads to pronounced gradients hampering the evaluation of the experimental data and even to irreversible decomposition of POM catalysts.

Conclusions
The present contribution investigates the transfer from batch to continuous operation under uncertainties for the POM catalyzed aerobic oxidation of DHA as an example for gas-liquid reactions. The study is based on a mathematical model assuming film theory, which consists of dimensionless balance equations for the liquid film and the liquid bulk at meso and macro scale, respectively. The parameterization is performed using experimental results obtained in a bSTR during the early stage of catalyst development. Importantly, the discrimination of the bSTR and the cPFR is realized by means of characteristic dimensionless numbers using the identical mathematical model for comparability reasons. The uncertainties typical for an early stage in process development are accounted for by sensitivity analyses carried out for the reaction kinetics and the mass transport parameter T.
It was revealed that the maximum degree of catalyst reduction strongly correlates with the mass transport parameter, but not directly with the reactor type. This is expected from the identical set of balance equations used and the re-oxidation of the catalyst chosen to be the reference reaction. The mass transport parameter, though, can directly be linked to characteristic values for both types of reactor operation and can thus be used for respective discrimination. The results also show that the mass transport parameter is suitable to determine, whether the reactor is operated under (nearly) intrinsic conditions or in a regime, where mass transfer, catalyst re−/oxidation and chemical reaction are interfering. Importantly, even rather slow reactions may suffer from mass transfer limitations in a bSTR and therefore require cPFR experiments for kinetic studies, which can be concluded from our results.
For the POM catalyzed aerobic oxidation of DHA we identified the catalyst reduction degree as the most important aspect to be considered in the design of the process. It appears to determine the significance of the experimental results with respect to catalyst activity improvement obtained under nonintrinsic conditions. Therefore, catalyst development needs to aim at sufficiently fast re-oxidation of the reduced catalyst species in order to avoid limitations or even catalyst decomposition. Furthermore, oxygen availability in the liquid phase supports catalyst re-oxidation and can be improved by enhancing the interfacial surface area, as well as higher oxygen solubility by the choice of appropriate solvents. Finally, the results indicate that determining the catalyst oxidation state in-situ is of superior importance in addition to the analysis of the product spectrum, in order to quantify the catalyst performance.
The study illustrates that even simple reactor models with roughly estimated parameters from preliminary experiments are sufficient in the early phase of process development. Such models are capable to resolve the complex interactions between mass transfer and reaction steps and thereby allow to identify the bottle-neck in the significance of experimental data. They can also be used for design of appropriate experiments and even for exploration of the feasibility of novel process windows. Hence, process development profits from both the transfer of experimental studies from batch to continuous operation and mathematical modeling in an early phase, already.