Novel Multiplicity and Stability Criteria for Non-Isothermal Fixed-Bed Reactors

With the increasing need to utilize carbon dioxide, fixed-bed reactors for catalytic hydrogenation will become a decisive element for modern chemicals and energy carrier production. In this context, the resilience and flexibility to changing operating conditions become major objectives for the design and operation of real industrial-scale reactors. Therefore steady-state multiplicity and stability are essential measures, but so far, their quantification is primarily accessible for ideal reactor concepts with zero or infinite back-mixing. Based on a continuous stirred tank reactor cascade modeling approach, this work derives novel criteria for stability, multiplicity, and uniqueness applicable to real reactors with finite back-mixing. Furthermore, the connection to other reactor features such as runaway and parametric sensitivity is demonstrated and exemplified for CO2 methanation under realistic conditions. The new criteria indicate that thermo-kinetic multiplicities induced by back-mixing remain relevant even for high Bodenstein numbers. In consequence, generally accepted back-mixing criteria (e.g., Mears’ criterion) appear insufficient for real non-isothermal reactors. The criteria derived in this work are applicable to any exothermic reaction and reactors at any scale. Ignoring uniqueness and multiplicity would disregard a broad operating range and thus a substantial potential for reactor resilience and flexibility.


INTRODUCTION
Currently, we see many incentives for more sustainable chemicals and energy carrier production based on CO 2 and H 2 . Chemical reactors for CO 2 hydrogenation play a crucial role in setting up sustainable production chains (e.g., via Fischer-Tropsch synthesis, CO 2 methanation, CO methanation, methanol synthesis, reverse water-gas shift). Especially for chemical energy conversion systems, as currently evaluated in the context of Power-to-X, these reactors are decisive for the overall process efficiency. Moreover, sustainable production increasingly demands for more flexible usage of chemical reactors accessing broader operation ranges and faster load changes. Encountering these new challenges demands a reassessment of former perspectives on reactor design and operation.
Due to the exothermic nature of most CO 2 hydrogenations, strong heat releases may result in harmful temperature excursions within the reactor interior, influencing catalyst lifetime, process safety, and performance. Cooled fixed-bed reactors allow for effective heat management and better controllability (Bremer and Sundmacher, 2019). These polytropic reactor concepts are often the first choice for the hydrogenation of CO 2 (Kopyscinski et al., 2010;Wood et al., 2012;Biollaz and Schildhauer, 2016).
However, designing and operating polytropic fixed-bed reactors is a rather challenging task, due to a variety of possible physiochemical phenomena (e.g., runaway, parametric sensitivity), which can lead to performance losses or critical process failures (see e.g., Kreitz et al., 2019;Kiewidt and Thöming, 2019;Moioli et al., 2019;Theurich et al., 2019;Fache et al., 2020;Fischer and Freund, 2020;Zimmermann et al., 2020). Such instances correlate with features from systems theory, i.e., uniqueness, multiplicity, and stability. Much work has been done in this direction for ideal plug flow tubular reactors (PFTRs) and continuous stirred tank reactors (CSTRs), providing explicit criteria based on fundamental reactor parameters (Zaldívar et al., 2003;Szeifert et al., 2007;Dochain, 2018). However, real polytropic reactors with finite back-mixing are yet inaccessible for these criteria and are therefore usually subject to qualitative descriptions.
This work intends to review the current perspective on uniqueness and multiplicity of real polytropic reactors under non-isothermal conditions and demonstrates their close connection to stability. Therefore, we perform a bifurcation analysis on a cascade of CSTRs with unlimited stage numbers. As a result, novel criteria for uniqueness and multiplicity are derived that underline the importance of back-mixing within fixed-bed reactors, even at an industrial scale. The theoretical discussions are exemplified for CO 2 methanationa highly exothermic reaction that is currently under intensive investigation for future energy storage.

Carbon Dioxide Methanation
The key reaction for methanation of CO 2 and H 2 reads Methanation is one of several possibilities to activate the inert carbon dioxide. Further reactions are, for instance, dry reforming, reverse water-gas shift (RWGS), and methanol synthesis. Among these reactions, CO 2 methanation is thermodynamically the most favored, as it exhibits the lowest Gibbs free energy of reaction in a temperature range from ambient to 900 K (Δ RG 0 −142 kJmol −1 ) (Wenzel, 2018). As depicted in Reaction 1, CO 2 methanation also features a strong exothermicity, such that lower temperatures shift the chemical equilibrium to the product side. This fact is supported by Figure 1, which shows the result of chemical equilibrium calculations assuming stoichiometric feed over a wide temperature range and for technical relevant pressures (thermochemistry data taken from Lemmon et al., 1998;Haynes, 2017). Figure 1 also illustrates two relevant kinetic models for nickel catalysts. The model from Koschany et al. (2016) differs significantly at temperatures above 700 K and shows higher methane contents due to the exclusion of carbon monoxide. Xu and Froment (1989), in contrast, considered the CO methanation and RWGS reaction to account for carbon monoxide. This was certainly motivated by the lower activity, which achieves sufficient conversions only at higher temperatures where CO formation becomes relevant. Although CO methanation and RWGS reaction certainly occur, they are not favored at elevated pressures and temperatures below 800 K. An analysis of thermodynamic limitations is comprehensively illustrated by Gao et al. (2012), who also showed possible amounts of solid carbon under CO 2 excess. More detailed calculations are presented by Kiewidt (2017), who identified considerable amounts of solid carbon at feed ratios H 2 /CO 2 < 3.2.
Regardless of the many details on chemistry, this work considers methanation mainly as a well-studied and highly relevant reference reaction. However, the presented results and criteria are easily applicable to other exothermic reactions.

Reactor Concepts for Exothermic Reactions
Reactor concepts for exothermic reactions are mainly determined by the expected temperature rise and the corresponding heat generation. In order to assess the temperature increase, the adiabatic temperature rise (ATR) is considered as an appropriate worst-case estimate. The ATR is defined as and exemplary evaluated for methanation at standard, stoichiometric feed conditions ΔT 0 ad 0.845 164'900 Jmol −1 0.044 kgmol − 1 2'950 Jkg − 1 K − 1 1073.5 K.
Looking at Eq. 2 reveals that the ATR is pressure-independent, since Δ RH and c p,gas depend exclusively on temperature under ideal gas conditions. Although a higher pressure leads to more reactants and thus more heat generation, the gas heat capacity increases simultaneously and compensates for a further temperature rise. The temperature dependence of the ATR is also rather negligible.
Technologically, this large temperature increase is answered with a series of adiabatic or polytropic fixed-bed reactors (with or without product recycles). Figure 2 illustrates exemplary state-space diagrams for both concepts. Which concept is used for a particular application is often driven by several aspects (e.g., plant size, product quality, costs, safety, reliability, flexibility). However, for modern hydrogenation, we see an increasing application of polytropic fixed-bed reactors (Kopyscinski et al. 2010).
Switching the reactor concept from adiabatic to polytropic (e.g., by cooled multi-tubular bundle reactors) potentially enables a lower catalyst volume, fewer reactor stages, better heat management, and the elimination of gas recycling. Nevertheless, such reactors are more difficult to construct and maintain (e.g., in case of catalyst replacements), which is, in particular, difficult for process scale-up. Hence, the performance benefit must outweigh the increased complexity of the reactor. In the context of gas-to-liquids, Fischer-Tropsch synthesis is an example of the successful use of multi-tubular reactors on a large scale (Wood et al., 2012). Figure 2 also highlights the relevance of uniqueness. Under certain conditions, both adiabatic and polytropic reactor concepts allow for nonunique operating states (here illustrated for CO 2 methanation with respect to cooling or inlet temperature). Especially for polytropic reactors, these nonunique states allow for large conversions and reduced peak temperatures in one single reactor stage, even for undiluted feeds (Bremer and Sundmacher, 2019). Furthermore, uniqueness is very much related to potential reactor runaways and is, thus, relevant to guarantee safe operation. Consequently, a clear determination of uniqueness is of great value for the design and operation of real non-isothermal fixed-bed reactors but has not been sufficiently addressed in the literature yet. This work shall contribute to close this gap and proceeds with a brief literature overview on statespace multiplicity.

State-Space Multiplicity
State-space multiplicity of fixed-bed reactors is differentiated in extrinsic and intrinsic (Bremer and Sundmacher, 2019), whereas this work exclusively deals with the intrinsic version. There are numerous sources for intrinsic state-space multiplicity. Three scales are in particular of interest: the scale of the active site, the catalyst pellet (including pores), and the reactor (Sheintuch, 1987). Figure 3 illustrates all three scales and the corresponding multiplicity sources (A to E) that have been most discussed in the literature.
These sources typically relate to mass and energy fluxes within a reactor control volume (zV) and their natural drive to equalize temperature and concentration differences caused by reactive sources and sinks. All fluxes are determined by corresponding transport resistances and linked in a network connecting scales and phases. Depending on inlet and surrounding conditions of the control volume, some resistances are rate-determining for the overall mass and energy transport. According to conservation laws, at steady state all fluxes, sources, and sinks are in equilibrium. A change of inlet and surrounding conditions disturbs the entire network, and a new flux equilibrium emerges. At reactor scale, there exists a unique flux equilibrium if all individual fluxes are unique and if inlet and surrounding of the control volume are not influenced by the control volume itself. If the control volume influences its inlet and/or surrounding (e.g., due to back-mixing), multiple flux equilibria, and, thus, multiple steady states are possible. This interaction is typically denoted as feedback at reactor scale. Similarly, these feedbacks may also occur at catalyst and site scale due to the interaction with neighboring scales.
When steady-state multiplicities in fixed-bed reactors are observed (either numerically or experimentally), they often result from several sources simultaneously. In order to assign the observed multiplicity to the correct source, careful evaluation of each source is required. Therefore, Sheintuch (1987) divided intrinsic multiplicity sources from Figure 3 into three classes: • purely kinetic (C m 2 ), • isothermal (A m , C m 1 ), • thermo-kinetic (A e , B e , C e 1 , D e , E e ).
Purely kinetic sources belong to mass transport at the active site induced by nonlinear kinetics (e.g., adsorption vs. reaction) (Nibbelke et al., 1998), concurrent reactions/reaction networks (Balakotaiah and Luss, 1983;Elnashaie and Elshishini, 1993;Nibbelke et al., 1998;Mohl et al., 2001), or catalyst deactivation (Gilles, 1977;Eigenberger, 1983). Isothermal sources arise from mass transport at catalyst (diffusion) (Morbidelli et al., 1986;Lee et al., 1987) or reactor scale (convection, dispersion). Thermo-kinetic sources are mainly caused by nonisothermicities, which influence the nonlinear reaction rate coefficients according to the Arrhenius relation (Eigenberger, 1972a;Heinemann and Poore, 1982;Wagialla and Elnashaie, 1995). Furthermore, multiplicities are also able to propagate through different scales and along with spatial directions. As a result, a large amount of multiple steady states (in theory up to infinity) at reactor scale can be produced by only a few (typically three) multiple steady states at pallet or site scale (Lee et al., 1987;Arce and Ramkrishna, 1991;Nibbelke et al., 1998). Due to the variety of different sources and their ability to propagate through scales, literature (going back to the early 60s) reports very differently or even controversially about total number, relevant sources, and necessary conditions of multiple steady states in fixed-bed reactors. This aspect is demonstrated by a brief literature review in Table 1.
As shown in Table 1, the investigation of steady-state multiplicity is mainly performed via cell models and dispersion models (see Materials and Methods). Both model types have been used to outline the two theoretical limiting cases: a reactor without and with an infinite number of multiple steady states. The dominating opinion is that the ideal PFTR is free of any state-space multiplicity due to the absence of transport in countercurrent direction (feedback). However, some studies opposed that the ideal PFTR is rather characterized by an infinite number of steady states (Liu and Amundson, 1962;Wagialla and Elnashaie, 1995;Nibbelke et al., 1998). These studies refer to the fact that the continuum description of an ideal PFTR corresponds to a series of infinite CSTRs. Assuming that the ignition can occur at any CSTR within the series gives rise to an ignition possible at any position along the reactor axis. This thought experiment allows for discontinuous solutions in packed-beds, which is often correlated to multiple steady states of single particles (Arce and Ramkrishna, 1991) or active sites (Nibbelke et al., 1998). Although the confusing concept of infinite solutions was discussed very controversially, it was certainly the main FIGURE 3 | Flux network and sources A to E for intrinsic multiplicity within fixed-bed reactors, double arrows indicate mass (m) and/or energy (e) fluxes to be part of the respective source.
Frontiers in Energy Research | www.frontiersin.org January 2021 | Volume 8 | Article 549298 driver behind many investigations. Thereby, two objectives became the most relevant for technical applications: Firstly, the ability to operate at various states promises performance improvements. Secondly, state-space uniqueness is of great value for safety reasons (e.g., to avoid runaways). Apart from ideal reactors, many studies in Table 1 also focus on real reactors considering finite Bodenstein numbers (see definition in Materials and Methods) within dispersion models or finite numbers n of representative CSTRs within cell models. As shown in Levenspiel (1999), both concepts are interchangeable due to the relation Besides the different model concepts, various assumptions are made in order to uncover the influence of specific sources of steady-state multiplicity. Therefore, the studies in Table 1 differ in many aspects of the model constitution (e.g., heterogeneous vs. homogeneous, isothermal vs. non-isothermal, first-order reaction vs. Langmuir-Hinshelwood Hougen-Watson-type reaction, with vs. without radial dispersion). For instance, the use of isothermal reactor models eliminates all thermo-kinetic multiplicities and allows for investigations on purely kinetic and/or isothermal multiplicities. On the contrary, the study of purely thermokinetic multiplicities is preferably done in the absence of kinetic and isothermal multiplicities.
Most studies dealing with thermo-kinetic multiplicity sources agree that axial dispersion of heat plays a vital role in the existence of multiple steady states in fixed-bed reactors. Eigenberger (1972b), in particular, argued that heat conduction through the solid phase (source B e ) must reduce an infinite multiplicity to a finite number of multiple steady states. The author also identified a maximum number of three or five multiple steady states, depending on the boundary conditions of the solid phase. It took about 20 years before this finding was revised. Therefore, Nibbelke et al. (1998) extended the model of Eigenberger (1972b) and considered a reaction kinetic with multiplicities at the active site (source C m 2 ). The authors could prove that an infinite number of steady states is maintained even if axial heat dispersion through the solid phase is incorporated. Unfortunately, this research direction is still very narrow, which is certainly due to the fact that the required heterogeneous model is a rather sophisticated  (2010) Purely experimental studies-comprehensive reviews given by Padberg and Wicke (1967), Wedel and Luss (1984), Harold and Luss (1985), Adaje and Sheintuch (1990) 2 exp (stable) Puszyński and Hlavacek (1984), Adaje and Sheintuch (1990) 4 exp (stable) Harold and Luss (1985) m -mass; e -energy; g -gas; c -catalyst; exp -experimental; k -kinetic; i -isothermal; tk -thermo-kinetic; CM -cell model; DM -dispersion model; MSS -multiple steady states.
Frontiers in Energy Research | www.frontiersin.org January 2021 | Volume 8 | Article 549298 numerical tool for such analyses. The use of pseudohomogeneous models provides a possible remedy. Although pseudo-homogeneous models are not able to distinguish between energy transport in source A and B, their results point in the same direction. By making use of dimensionless model formulations and bifurcation techniques Heinemann and Poore (1981), Jensen and Ray (1982), Heinemann and Poore (1982) classified several aspects that also contribute to the existence of multiplicity. Both stated that a nonadiabatic reactor with unequal Bodenstein numbers for heat and mass dispersion shows a unique solution, either for sufficiently high values of the Bodenstein numbers, large heat transfer coefficients, or small values of the Damköhler number. These conditions also hold for industrial fixed-bed reactors; such that multiplicity was found to be relevant even for technical configurations (Puszyński et al., 1981;Pita et al., 1989). This is reflected in Table 1 by the broad range of Bodenstein numbers in which multiple steady states are observed. In most investigations, the number of three multiple steady states is confirmed. Beyond that, some studies report up to eleven multiple steady states. However, experimental evidence above four multiple steady states was not found. For instance, Wedel and Luss (1984) validated their results with an experimental setup for CO and CO 2 methanation in a fixed-bed of 25 cm in length and 2.5 cm in diameter. The authors could reproduce predicted model solutions of one ignited and one extinguished state with identical feed concentration. Therefore, a one-dimensional axial dispersion model, which accounts for the diffusion-reaction interactions within the pellets, was found to be well suited for multiplicity analysis. Since their experimental setup has a low bed length and, thus, a low Bodenstein number (Bo e 8), conclusions for industrial reactors are rather limited. Nevertheless, the good agreement between model and experiment verified and proved that state-space multiplicity is of significance within fixed-bed reactors. More recently, Agrawal et al. (2007) pointed out that the often reported high number of multiple solutions originates from the assumption of constant transport coefficients, which leads to many fragile solutions that emerge from unstable branches. Accounting for the variation of heat and mass transfer coefficients with local flow and reaction properties eliminates these nonphysical solutions. For this reason, recent works are often motivated to add more physical meaning to their models. This allows to additionally focus on effects of secondary order (e.g., flow maldistributions, localized hot-spots, spatial and spatiotemporal patterns) (Sheintuch, 1997;Trinh and Ramkrishna, 1997;Jaree et al., 2001;Papadias et al., 2001;Marwaha and Luss, 2003;Agrawal et al., 2007;Viswanathan et al., 2008;Nekhamkina and Sheintuch, 2012), which is, however, not object of this work. Similarly, stable oscillatory solutions are also disregarded in this work, since they are unlikely to occur in fixed-bed reactors on an industrial scale due to their high thermal inertia (Jensen and Ray, 1982).
In order to identify the full operating range of fixed-bed reactors, this work considers state-space multiplicity as a crucial feature. Therefore, thermo-kinetic sources are of major interest, since purely kinetic and isothermal multiplicities are reported to be rather fragile. The literature survey from above reveals that axial dispersion, as well as the diffusion-reaction interaction, are essential model components to obtain consistent results. However, for specific applications, these studies rather provide trends than generally valid correlations and criteria. The impact of several simultaneously interacting features is certainly the main reason for this obstacle. Furthermore, the operation at unstable states is also underrepresented in the literature and needs further elaboration (see e.g., Bremer and Sundmacher, 2019). Therefore, the following sections intend to provide more insights and propose general criteria that unifies the aspects uniqueness, multiplicity, and stability.

MATERIALS AND METHODS
Analyzing chemical fixed-bed reactors requires to incorporate a variety of physicochemical interactions and is often motivated by models from first-principles. Solving these models can lead to a challenging computational task due to the complex coupling of mass, energy, and momentum transport. However, putting all available details into consideration is not always necessary. For instance, the basic phenomena described in this work are fully representable by mass and energy transport alone, which agrees well with previous works (see Table 1). In addition to physical assumptions, the considered control volume boundary also determines the relevant balance components. This differentiation meets the two ideal reactor concepts for continuously operated reactors: • the ideal CSTR, • the ideal PFTR.
The CSTR concept often considers a control volume around the entire reactor volume (integral balance). In contrast, the control volume for the ideal PFTR typically refers to an infinitesimally small volume segment within the reactor volume (differential balance). The PFTR allows to describe in detail spatial distributions of the reactors state variables (e.g., temperature and mass fraction), but also requires higher numerical efforts. Both ideal reactors also represent the two limiting cases for back-mixing/axial dispersion, which is typically measured by the axial Bodenstein number for mass and energy according to which are essential for later discussions. In case of an ideal CSTR the reactive volume is fully mixed (D z λ z → ∞, Bo m Bo e → 0), whereas in case of the ideal PFTR no backmixing/dispersion exists (D z λ z → 0, Bo m Bo e → ∞) (Levenspiel, 1999).
Real reactors, as considered in this work, are allocated right in between these two limiting cases. The corresponding model approaches are • the CSTR cascade (tanks-in-series model, cell model), • the tubular reactor model with axial dispersion (dispersion model, continuous model).
In order to represent real reactors more accurately, the control volume might also differentiate between bulk gas phase and catalyst phase. Two model concepts are commonly applied: • the pseudo-homogeneous model (no phase distinction), • the heterogeneous model (phase distinction with interfacial and intraparticle mass and energy transport).
Depending on the considered model, the energy dispersion coefficient in Eq. 4 is considered differently. The pseudohomogeneous model lumps solid and gas phase, which requires an effective energy dispersion term (λ z λ eff z ). On the contrary, the heterogeneous model allows for the distinction between gas and catalyst-phase dispersion and, thus, accounts for two separate energy-based Bodenstein numbers (Bo eg , Bo ec ). For real fixed-bed reactors, the mass-based Bodenstein number is reported to be three to ten times higher than the energy-based Bodenstein number (Puszyński et al., 1981;Bostandzhiyan and Shkadinskii, 2010). Thus, energy back-mixing is the dominating axial dispersion mechanism within fixed-bed reactors (Mears, 1976). Note that some studies define the here stated Bodenstein number as Péclet number. In this work, these numbers are distinguished by the characteristic length, which is the reactor length L for the Bodenstein number and the particle diameter d p for the Péclet number.
Both model concepts offer different levels of sophistication, e.g., one-, two-, or three-dimensional spatial resolution, stationary or dynamic, with or without axial dispersion. How sophisticated a model needs to be is often rated by criteria (e.g., Mears' criterion), which incorporate dominating transport phenomena and reactor design. A comprehensive overview is provided by Pérez-Ramírez (2000). Besides these criteria, distinct model or transport components may be examined separately. This decoupling approach is well-established to identify and analyze, for instance, different sources of state-space multiplicity. In this context, Nibbelke et al. (1998) emphasized that the identification of multiplicity sources is inevitable for a correct interpretation of numerical and also experimental results. This observation is crucial for the construction of adequate reactor models as well as for optimal reactor operation and control.
In the following, a single pseudo-homogeneous CSTR model is used to represent fully mixed reactive conditions and to illustrate its implications on uniqueness, multiplicity and stability. Afterward, the single CSTR is extended to a series of CSTRs (cell model) in order to approach conditions as they prevail in real fixed-bed reactors. All investigations are exemplified for methanation under realistic reactive conditions.

RESULTS
As indicated by the previous section, back-mixing is an essential feature that determines multiplicity. In order to show how this applies to real reactors, the following derivations shall serve as a guideline for multiplicity analysis, exemplified for catalytic methanation. Beginning with the well-known limiting case of a fully mixed reactive volume illustrates the fundamental parameters that characterize the range and intensity of multiplicity.

CSTR Analogy
The technical relevance of this limiting case can be found, for instance, with Berty reactors, which are widely used for catalyst characterization. Although the relationships presented here are already state of knowledge, they are crucial for later excursions toward industrial-scale fixed-bed reactors. Further on, a CSTR model as described in Supplementary Section S1 is used. The required mass and energy balances read accordingly: Consequently, steady-state operation of a single CSTR is governed by the equality of Eqs. 5 and 6. In order to evaluate unique and non-unique operating points, Supplementary Table S1 provides a reference setting that corresponds to carbon dioxide methanation in a fixed-bed reactor including realistic parameter ranges.
The six key parameters in Supplementary Table S1 are highlighted in bold and result in Figure 4 for variations of the operating temperature and Figure 5 for variations of the Stanton number. The operating temperature is given in Supplementary Eq. S19 and reflects changes in coolant and inlet temperature, whereas the Stanton number mainly results from changes in the coolant heat transfer. The left side of both figures shows that under realistic conditions, multiple equilibrium (or operating) points are indeed attainable, similar to the theoretical discussions in State-Space Multiplicity. The right side in both figures illustrates how hysteresis emerges when the operating temperature varies within realistic ranges. Apart from variations in operating temperature and Stanton number, one might also consider variations of the residence time τ. Higher residence times correspond to higher Damköhler numbers and increase the curvature of the mass balance operating curve (Eq. 5) and, thus, increase the hysteresis.
These results explicitly show that under perfect back-mixing at most three operating points are attainable, whereas two are stable (OP 1 and OP 3) and one is unstable (OP 2). As seen in the upper right figures, the unstable operating points are always surrounded by stable ones and cover a significant part of the attainable conversion range. If these unstable states are ignored for reactor design and operation, a significant potential might get lost. This clearly demonstrates the close connection between multiplicity and stability. Here, the generalized criterion for stability of a CSTR at a certain operating point is known to result from mass and energy-based sensitivities according to which is equivalent to the criteria proposed by van Heerden (1953) and Gilles and Hofmann (1961). Note that sensitivities might also be considered with respect to other parameters than temperature. However, this work focuses primarily on thermo-kinetic multiplicities, for which temperature sensitivities are the most relevant. The mass-based sensitivity may be further evaluated by resolving the total differential at constant τ and c in , such that Plugging this into Eq. 7 leads to the rate-based CSTR stability criterion CSTR stability : In comparision, Szeifert et al. (2007) comprehensibly reviewed available criteria for reactor stability, runaway, and parametric sensitivity. Among those, the authors recommend using the Ljapunov stability criterion for a PFTR in space (or for a batch reactor in time), which is very similar to the upper CSTR stability criterion. In case of an equimolar reaction the PFTR stability criterion accordingly reads PFTR stability : Both criteria reveal that reactor stability is supported under intensive cooling conditions (St → ∞) and suppressed under reduced cooling conditions (St → 0). As explained in State-Space Multiplicity the PFTR is primarily considered to have a unique solution for fixed parameters. Hence, the features stability, runaway, and parametric sensitivity need to be separated from multiplicity, which was already highlighted by (Bilous and Amundson, 1956). Until today, this separation is often considered to be generally valid. However, in the CSTR case a strong coupling of these features is clearly given. From the above explanations, generalized criteria for uniqueness and multiplicity within the operating range O derive very similar to the stability criterion and read CSTR uniqueness : max CSTR multiplicity : max On the left side, the mass-based sensitivity corresponds to the steepest conversion gradient at the inflection point T T * of Eq. 5. On the right side, the energy-based sensitivity corresponds to the constant conversion gradient of Eq. 6. Although the uniqueness criterion applies to all operating points in O (sufficient condition for uniqueness), the multiplicity criterion only indicates the existence of some operating points with multiple steady states in O (necessary condition for multiplicity). Thus, even if the multiplicity criterion holds, unique solutions in O might still exist. Furthermore, discontinuous jumps at the turning points from a stable extinguished to a stable ignited branch -and vice versa -is also solely induced by the existence of multiplicity. These jumps coincide with what is generally denoted as parametric sensitivity and reactor runaway (Morbidelli and Varma, 1982). Moreover, the equivalence strikingly demonstrates the close connection between stability, runaway, parametric sensitivity, and multiplicity when backmixing is present. However, commonly accepted criteria for reactor runaway and stability do not consider back-mixing at all (Szeifert et al., 2007;Kummer and Varga, 2019). To close this gap, the next section applies the upper correlations for mass and energy to a series of CSTRs, where a finite number of stages corresponds to finite back-mixing.

CSTR Cascade Analogy
In order to approach a more realistic description of industrialscale fixed-bed reactors, this section extends the previous relations to a CSTR cascade with multiple stages. Note that such cascades also mimic finite volume schemes, which are widely used to solve dispersion models numerically.
Supplementary Section S2 provides the entire cascade model notation, mass and energy balances, as well as the solution strategy. The resulting equations of each CSTR stage {i} are In order to separate different multiplicity sources, isothermal and non-isothermal cascades are discussed in the following.

Two Isothermal CSTRs
Under intensive cooling conditions (St → ∞) or negligible heat effects (ΔT {i} ad → 0) the CSTR cascade becomes isothermal, meaning that T in T cool T {i} . The corresponding implicit nonlinear equation system results from Eq. 11 and reads which can be solved consecutively (stage-wise) or simultaneously for X {1} CO2 and X {2} CO2 by root-finding algorithms (e.g., Newton-Raphson method). The solution of Eq. 13 is illustrated in Figure 6.
As seen on the left side of Figure 6, the operating point (depicted by bullets) of each CSTR lies on a straight line in the R-X diagram. The unique intersection at any temperature proves, that no purely kinetic multiplicity (C m 2 in Figure 3) exists for the rate expressions used in this work. Since intraparticle and interfacial transport resistances are neglected (η meth 1), the multiplicity sources C m 1 and C e 1 are also absent. This graphical analysis is also applicable to other reactions in order to check for these multiplicity sources. The right side of Figure 6 depicts the stage operating range as well as the overall operating range for two CSTRs in series and compares it to scenarios with 1, 10, and 100 stages. The comparison shows the well-known tendency that the higher the stage number, the higher the conversion.

Two Non-Isothermal CSTRs
The operating points for a series of non-isothermal CSTRs are additionally determined by individual stage temperatures deviating from inlet and cooling temperature T in T cool ≠ T {i} . Thus, the nonlinear equation System 13 needs to be extended by the stage energy balances in Eq. 12 such that This system can be solved again stage-wise or simultaneously for X {1} CO 2 , X {2} CO 2 , T {1} , and T {2} by root-finding algorithms. The previous section showed that one single stage allows for up to three multiple solutions, which indicates that the CSTR cascade in Eq. 14 also generates multiple steady states. In theory, if each stage exhibits up to three multiple steady states (OP 1, 2, and 3), a total number of 3 n state combinations arise for the entire cascade. Standard root-finding algorithms, however, converge only to one local solution of System 14. Thus, finding all solutions requires further efforts, which is part of the following discussions.
The graphical solution of System 14 is shown in Figure 7 considering three of nine possible combinations ([1-1], [2-1], ). This figure comprehensively shows how the first operating point influences its downstream stage. The more conversion is achieved within the first stage, the less remains within the next stage. Furthermore, the energy-based operating line (green) of the second stage becomes steeper due to the reduced ATR. Both effects simultaneously reduce the occurrence for multiplicities within the second stage, if the first stage operates at an ignited state (OP 2 or 3). Later, this will be an essential aspect to interpret solutions of multi-stage CSTRs. Enumerating all possible solutions within the relevant coolant temperature range leads to the results in Figure 8.
Compared to the single CSTR, the hysteresis is more pronounced, and intermediate solutions emerge. Although combinatorics allows nine multiple steady states, only seven remain significant after enumeration. The combinations [3-2] FIGURE 6 | State-space diagram for a cascade of two isothermal Continuous Stirred Tank Reactors (CSTRs) in the R-X (left) and X-T (right) plane, reference setting taken from Supplementary Table S1 but T op 550 K.  and [3-3] cancel out, due to the previous discussions on Figure 7. In addition, combinations [2-2] and [2-3] are rather seen as fragile because they exist in a very narrow operating range. In consequence, multiplicity driven by combinatorics is divided into significant and non-existing paths, as illustrated in Figure 9. These findings indicate that multiplicity in a CSTR cascade attributes only to a few stages. In order to confirm this hypothesis, a multi-stage cascade is analyzed in the following.

Multi-Stage CSTR Cascade
The enumeration effort of a multi-stage CSTR cascade increases exponentially with the number of stages if all possible combinations are taken into account. With six stages, there exist already 3 6 729 state combinations for all coolant temperature increments (here 300), leading to over 200'000 solver runs. Although a brute force enumeration would provide all solutions, many of them will be insignificant, due to missing physical constraints. Thus, a more elegant way is the use of bifurcation theory and numerical continuation techniques, which track solution branches according to slight parameter changes (e.g., coolant temperature) within a predefined range. During these changes, the eigenvalues of the linearized system equations may change such that the system becomes unstable. The shift to instability occurs if some eigenvalues cross the imaginary axis. At this point, bifurcation can take place from which new branches spread. However, this method also suffers from extensive computational efforts due to the eigenvalue calculation and possible inaccuracies of the systems Jacobian close to the bifurcation point. Furthermore, there is no guarantee that this method obtains all solutions. More details on bifurcation theory and numerical continuation techniques can be found in (Heinemann and Poore, 1981;Jensen and Ray, 1982;Kubíček and Marek, 1983;Wagialla and Elnashaie, 1995).
In this work, the enumeration of all solutions was found to be the more convenient and illustrative approach. It will be shown that the non-existent solutions can be excluded already in advance. Together with efficient nonlinear computation techniques (here provided by CasADi; Andersson et al., 2018) the enumeration of all solutions with reasonable computational effort becomes feasible. Accordingly, the results for three to six CSTRs in series are illustrated in Figure 10.
Most importantly, the solutions in Figure 10 aggregate in three main clusters forming -similar to a single CSTR -a stable ignition and extinction branch, as well as an unstable intermediate branch. The number of multiple steady states increases with increasing stage number n, but most solutions persistently converge to the three main clusters. Once again, each relevant solution has not more than two stages exhibiting multiplicity. Cascades with more stages as represented by Figure 11 confirm that two stages are sufficient to map the dominating state clusters. These two stages are further on denoted as key stages, which may exist at any position within the cascade.
The key stage solutions directly at the inlet are colored in Figure 10. Solutions with the same key stage combination but different locations along the cascade aggregate within the same cluster. For instance, the upper ignition branch is covered by n [/1-3-1/] combinations, the middle unstable branch is FIGURE 8 | Enumerative solution of significant operating points in a two-stage CSTR series, reference setting taken from Supplementary Table S1 but T in T cool . Note: T cool is constant over the entire cascade, whereas T op changes in each CSTR stage.  Figure 12, showing the relevant state combinations of a multi-stage CSTR cascade. As illustrated in Figure 11, the clustering of all key stage solutions still remains for higher stage numbers. It was found that the hysteresis loop widens significantly with increasing stage number. However, after reaching a certain stage number, the hysteresis loop degenerates again until it finally disappears. Furthermore, different shapes of the unstable intermediate branch emerge, depending on the parameter setting. The unstable branch evolves more or less pronounced, as indicated by variations of the residence time in Figure 11. This is in particular relevant if an operation at these unstable branches is aspired (Bremer and Sundmacher, 2019).
Technically, key stages are the analogy of narrow reaction fronts in real fixed-bed reactors. Those fronts often develop within a very short reactor segment, preferably close to the reactor inlet. Considering this, the key stages become less likely if they are located further downstream, which curtails the relevance of the black subordinate branches associated with the last stage. In real reactors, the actual reaction front position is typically determined by second-order effects (e.g., preheating, dispersion, flow maldistribution, heat conduction within the reactor jacket).
From all this, the following three-level hierarchy of thermokinetic multiplicity in fixed-bed reactors can be drawn:

Assumption
Multiplicity in all stages Two key stages reaction front Cluster formation Max #MSS 3 n 0 4n − 1 0 3 to 5 FIGURE 10 | Enumerative solution of operating points (OPs) in a multi-stage CSTR cascade, color-OPs of two key stages at the inlet, black-OPs of subordinate branches, reference setting taken from Supplementary Table S1 but T in T cool . Moving from left to right considers more physical details and approaches the situation in real fixed-bed reactors. This threelevel hierarchy compromises many different and controversial opinions found in the literature (see State-Space Multiplicity). For instance, the findings of Eigenberger (1972a), andEigenberger (1972b) reporting only three to five multiple steady states are mainly associated with cluster formation, whereas studies that report an infinite number of multiple steady states (for n → ∞) neglect clustering and count each state separately.
So far, state clustering and its connection to reaction fronts and multiplicity in real reactors is rather disregarded in literature. One reason might be the missing availability for efficient numerical tools, which have only been accessible in recent years. Since the previous considerations mainly address multiplicity trends, it remains to be shown under which conditions uniqueness applies.
Until today, an exact uniqueness criterion for nonisothermal fixed-bed reactors largely remains an open question (Dochain, 2018). The axial dispersion model has been a favorite target for extensive mathematical analyses but often limited to first-order reactions under isothermal conditions (Schmitz, 1975;Varma, 1980;Arce and Ramkrishna, 1991). However, numerous studies provide qualitative trends pointing in similar directions. For instance, Jensen and Ray (1982) summarized, that the solution will be unique for sufficiently high Bodenstein numbers, large heat transfer coefficients, or small Damköhler numbers. The previous results are very much in line with these qualitative trends, but moreover, they reveal generalized criteria for stability, uniqueness, and multiplicity of non-isothermal fixed-bed reactors. These criteria are derived in the following.

Stability, Uniqueness, and Multiplicity Criteria for Real Non-Isothermal Reactors
The observations from above show that the thermo-kinetic multiplicity feature of the first stage is entirely capable of representing the three main state clusters. Consequently, if the first stage is free of multiplicity, then all following stages are also free of multiplicity. This key feature enables the criteria of the first CSTR stage (see Eqs 8-10) to be assigned to the entire CSTR cascade (cell model) according to with T T * as inflection point of Eq. 11 corresponding to the steepest conversion gradient. For simplicity, the index CO 2 is omitted here and in the following.  Consequently, in the absence of purely kinetic and isothermal multiplicity, uniqueness in non-isothermal CSTR cascades boils down to very few key parameters, lumped together as mass and energy-based thermal sensitivity. Both sensitivities are evaluated in Figure 13 with respect to the methanation reference setting of Supplementary Table S1.
According to the reference setting in Figure 13, uniqueness can only be guaranteed for cascades with several thousand stages. This fact still applies to a wide range of heat transfer and catalyst activity, as indicated by variations of the Stanton number and effectiveness factor. Intensified heat transfer mainly affects cascades with higher back-mixing and leads to reduced multiplicity regions. In some scenarios, uniqueness is guaranteed for low and high back-mixing conditions, but not for the intermediates (e.g., for St 100 and η 1). In contrast, reducing the ATR (e.g., via product gas recycling) always leads to diminished multiplicity regions. At adiabatic conditions (St → 0), backmixing does not influence the energy-based thermal sensitivity, and uniqueness becomes solely determined by the mass-based thermal sensitivity and the ATR. Not shown is the influence of the remaining key parameters, pressure, and residence time. However, both are indirectly incorporated into the Stanton number and effectiveness factor.
The equivalence of stage and Bodenstein number in Eq. 3 at low back-mixing (Bo > 100) allows for the transition from cell to dispersion models. a Therefore, we introduce a surrogate conversionX which allows for applying the first stage mass balance from Eq. 11 to dispersion models according to which is then used to calculate the mass-based thermal sensitivity of a fixed-bed reactor. Together with the adapted energy-based thermal sensitivity from Eq. 12 at elevated Bodenstein numbers (Bo > 100) the previous criteria read as stability(DM) : with T T * as inflection point of Eq. 11 corresponding to the steepest conversion gradient. The evaluation of the reference setting from Supplementary Table S1 leads to the same results as in Figure 13 with n Bo/2. The generalized Criteria 16-18 have not been found in literature yet. They can be used as an a priori estimate for any exothermic reaction and reactors at any scale, only requiring apparent rate expression, coolant heat transfer coefficient, inlet condition, and back-mixing intensity. Note that no expensive computation of the entire dispersion model is required, which makes it easy to use for reactor design, operation, and safety analysis. The criteria, however, demand for a representative Bodenstein number either mass or energy-based. As illustrated in State-Space Multiplicity, the energy feedback is of major interest for thermo-kinetic multiplicity, which indicates that the three to ten times smaller energy-based Bodenstein number is the most reasonable choice. Furthermore, the here proposed criteria recommends considering axial dispersion even for high Bodenstein numbers beyond 400. This is contrary to commonly accepted criteria of Hlaváček and Marek (1966), Mears (1976), Young and Finlayson (1973), and Mederos et al. (2009), which did neither consider multiplicity nor reactor stability.
In order to access the mass-based sensitivity on the left, the total differential applied to the implicit Eq. 15 may help: Consequently, the sensitivity is represented by dX dT With Eq. 20, the mass-based sensitivity can be evaluated within the relevant temperature range to identify the maximum gradient at T * and Da * I . Thus, there is no further need to solve the implicit Eq. 15. Note that the here used Damköhler number remains as a function of temperature such that the above criteria condense to a This equivalence is exploited when the finite volume upwind scheme is used to solve dispersion models numerically. In this case, the corresponding number of finite volumes in flow direction inherently contains a certain degree of backmixing. If the number of finite volumes is too low, an artificial dispersion (so-called numerical diffusion) will superimpose other dispersion components included in the model. stability 1 st − order, equimolar : with T * and Da * I as arguments of the maximum gradient of Eq. 20. Although these simplified criteria are not adequately applicable for methanation due to the strong influence of the thermodynamic equilibrium, they comprehensively show how the key parameters affect multiplicity. At the beginning of this section, the current state in the literature was highlighted to be rather qualitative. The statement of Jensen and Ray (1982) saying that the solution will be unique for sufficiently high Bodenstein numbers (Bo → ∞), large heat transfer coefficients (St → ∞), or small Damköhler numbers (Da → 0) is perfectly represented by the Criteria 22 and 23. Furthermore, the limiting case of an adiabatic CSTR (St → 0, Bo/2 → 1) applied to Criteria 21 is equivalent to the stability criterion presented by Kimura and Levenspiel (1977).
In extension to the infinite back-mixing case in CSTR Analogy, it can be confirmed that stability and uniqueness are closely related under finite back-mixing conditions. However, the finite back-mixing case (1 < Bo/2 < ∞) is typically not considered in the literature (Szeifert et al., 2007;Kummer and Varga, 2019), or was found to be insignificant (Balakotaiah et al., 1995) for stability analysis. In contrast, this work shows a distinct relevance of multiplicity for real reactors, so that back-mixing must also be highly relevant for stability, runaway, and parametric sensitivity. In this regard, the quantitative description of all features results in the here proposed Criteria 17 and 22.

DISCUSSIONS
In summary, uniqueness and multiplicity of real non-isothermal reactors have proven to be decisive characteristics. In addition to the qualitative descriptions prevailing in the literature, the criteria proposed here represent a novel quantitative measure applicable to any exothermic reaction and reactors at any scale. This work also shows that mass and energy back-mixing represents the essential link between uniqueness, multiplicity, stability, runaway, and parametric sensitivity, which are usually treated independently in the literature. The observation that back-mixing remains relevant for these characteristics even at high Bodenstein numbers implies that generally accepted back-mixing criteria of Hlaváček and Marek (1966), Mears (1976), Young and Finlayson (1973), and Mederos et al. (2009) are insufficient for real non-isothermal reactors.
In addition, the CSTR cascade model indicates that a narrow reaction front, not larger than two representative CSTR stages (key stages), mainly determines uniqueness and multiplicity in real non-isothermal fixed-bed reactors. The illustrated stateclustering of the cascade model accounts for the fact that this reaction front may occur at any position within the fixed-bed. From this, a three-level hierarchy is derived, which unifies controversial opinions that still exist in the literature.
In summary, the presented methodology, as well as the derived criteria, shall allow for easier accessibility of fundamental reactor characteristics. This is particularly useful for the future objective of operating chemical reactors more flexible and within larger operating ranges.
FIGURE 13 | First stage mass and energy-based sensitivity for non-isothermal CSTR cascades with various stage number, effectiveness factor and Stanton number, reference setting taken from Supplementary Table S1.
Frontiers in Energy Research | www.frontiersin.org Therefore the derived criteria may be used for reactor design, control, and safety purposes. Ignoring uniqueness and multiplicity would disregard a broad operating range and thus a substantial reactor performance potential.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JB and KS contributed conception and design of the study. JB organized the content, ran the computational experiments, and derived the criteria. KS motivated the detailed bifurcation analysis and supervised the project. JB wrote the first draft of the manuscript. All authors contributed to manuscript revision, read and approved the submitted version.