The Formalism of Chemical Thermodynamics Applied to an Oscillatory Multistep Chemical System

The thermodynamic optimization of a process focuses on consumption, production, and efficient use of energy. The unsteady-state nature of batch reactor processing requires describing the set of processes’ dynamic behavior for energy optimization. This work aims to apply the formalism of chemical thermodynamics to a multistep chemical system in a batch reactor, aiming for a dynamic description of its evolution to the equilibrium state. As the system of study, we selected a mathematical model called the Oregonator , derived from the mechanism of the oscillating Belousov-Zhabotinsky reaction. In the methodology, we used the reaction quotient to evaluate the Gibbs function, the thermodynamic affinity, and the entropy generation as a function of the reaction extent. The results show that the overall reaction fulfills the thermodynamic fundamentals of chemical equilibrium, despite having a non-stoichiometric coefficient. However, the multistep coupled reaction system does not allow verifying compliance with the thermodynamic foundations of chemical equilibrium. We conclude that it is necessary to improve thermodynamic formalism to describe multistep chemical processes as a function of a global reaction extent variable. In this scenario, the entropy production rate emerges as a promising quantity.


Introduction
The study and characterization of chemical equilibrium are of great importance for describing processes that depend on the phase stability of fluids, such as separation in reactors and reactive and azeotropic distillation in columns (Prausnitz et al., 1998;Saunders and Miodownik, 1998;Tosun, 2021).On the other hand, in reacting systems, the chemical equilibrium defines the upper bound of transformation of a chemical reaction under the initial conditions and the imposed constraints.A closed system that evolves at constant T and p must satisfy the constraints of directionality and equilibrium associated with thermodynamic potentials (Rosenberg and Klotz, 2008;Honig, 2020).
Let us consider a reaction system composed of n simultaneous chemical reactions coupled to each other through intermediary chemical substances.The set of chemical equations that satisfies the law of mass conservation and describes the reactant system composed of m substances can be written compactly, as shown in Equation ( 1 The reaction extent r ξ , as a state variable of the reaction system, allows evaluating the change in the Gibbs function of a chemical system described by Equation (1) by following Equation ( 4): At constant T and p , and due to the initial conditions of the set of chemical reactions, the change in the reaction extent has an arbitrary sign.Therefore, the equilibrium condition , 0 T p dG = must be satisfied by the closed system of chemical reactions; based on Equation (4), this is described for each r reaction in Equation ( 5): ( ) The equilibrium condition in Equation ( 5) can also express a generalized force, i.e., the thermodynamic affinity A , according to Equation (6): A contracted form of Equation ( 6) can be written in terms of a reaction quotient and an equilibrium constant (Dutt, 1985).From Equation ( 4) and Equation ( 6), and due to the irreversible nature of chemical reactions, for a reacting mixture in a closed system at constant T and p the Gibbs function reaches a minimum at the chemical equilibrium defined by a thermodynamic constant.The evolution towards and at the equilibrium state of the system described by Equation (1) satisfies the conditions given in Equation ( 7).
Thus, as the reaction progresses toward the chemical equilibrium, the thermodynamic affinity decreases until it becomes equal to zero (Solaz-Portolés, 2011;Kondepudi and Prigogine, 2014).
The above-presented formalism allows for a direct account of the irreversible isothermal transformation of matter during a chemical reaction through the generation of entropy (thermodynamic dissipation) described in Equation ( 8) (Prigogine, 1961;Kondepudi and Prigogine, 2014;Barragán et al., 2021).
As per Equation ( 8), the thermodynamic dissipation is zero at the chemical equilibrium, since the affinity and the net reaction rate are zero.
Although the above is widely accepted and duly documented in many texts (Eu and Al-Ghoul, 2018), great doubts arise when applying it to concrete systems (Beretta and Gyftopoulos, 2015), such as those raised by the discovery of oscillating chemical reactions in the 1950s.By then, agreeing on chemical oscillations and the second law of thermodynamics was unclear (Kiprijanov, 2016).Two decades later, once experimentation validated the first observations of chemical oscillations and provided an understanding on how the second law of thermodynamics rationalizes these systems, a fruitful field of research began, known as nonlinear chemical dynamics (Epstein and Pojman, 1998).There is extensive development in the formalism and the applications of the chemical equilibrium of phases.However, there are still doubts about the chemical equilibrium in reacting systems, particularly its teaching and the conceptual errors that can affect it (Quílez-Pardo and Solaz-Portolés, 1995;Banerjee, 1995;Rogers et al., 2000;Patiño-Sierra and Barragán, 2022;Martínez-Grau et al., 2014).
In this work, we test the application of thermodynamics formalism to multistep reacting systems using the mathematical model known as the Oregonator, derived from the chemical mechanism of the Belousov-Zhabotinsky (BZ) reaction. (1) (2) (3) (4) (5) (6)

Methodology
This section describes the Oregonator model and the kinetic parameters used in this study.It presents the equations that govern the dynamics of the model's overall reaction and the system of equations for each of the reaction steps.We apply the formalism of thermodynamics to derive the expressions for the reaction coefficient and the Gibbs function.All equations are solved numerically by using an integration subroutine for systems of stiff differential equations, which is available in the MATLAB software.

Multistep chemical system: the Oregonator
The Oregonator model describes the main characteristics of the oscillatory dynamics of the BZ reaction, both in an open and closed system (Field et al., 1972;Field and Noyes, 1974).The Oregonator comprises five reaction steps, four for a set of inorganic reactions of the reaction BZ and one for organic reactions: Equations O1 to O5 in Table 1.The reaction Og is the overall reaction of the model, obtained as follows: Og=2O1+O2+O3+O4+O5.The overall reaction Og satisfies the following: A and B reactants have integer stoichiometric coefficients, intermediates X and Z are absent, and the stoichiometric coefficient f is uniquely related to Y .The variables in the Oregonator model correspond to the following chemical species in the BZ reaction: A to sodium bromate, B to malonic acid, Z to cerium (IV) ion catalysts, X to bromous acid, Y to bromide ion, and P to hypobromous acid.Equation O3 is the chemical instability of the model, which is the autocatalytic production of bromous acid.Equation O5 represents a set of organic reactions, with malonic acid, the organic substrate, kinetically controlling the concentration of molecular bromine and bromide ions, which are the species driving the positive feedback loop in the reaction.The stoichiometric coefficient f in equation O5 accounts for the complexity of the chemical reaction and acts as a bifurcation parameter in the Oregonator model (Field and Noyes, 1974).Although there are many versions of the Oregonator model, whose majority does not include the B variable, in this work, we employ the original version derived from the BZ reaction mechanism to use the organic substrate concentration as a dynamical control parameter (Table ) 1 (Tyson, 1981(Tyson, , 1982)).

Mathematical model: the overall reaction Og
When applying Equation (2) to the overall reaction Og, the change in each variable as a function of the reaction extent is as follows: The reaction rate, expressed as a function of the reaction extent by and applying the mass action law to Og, is shown in Equation ( 9).

Mathematical model: reaction steps O1 to O5
When applying Equation (3) to the system of equations O1-O5, the change in each of the variables, expressed as a function of the progress of each reaction r , is as follows: Thus, by applying the law of mass action, Equations ( 10)-( 14) show the reaction rate of the steps in the Oregonator model: Set of equations and kinetic constants of the Oregonator model (Tyson, 1981(Tyson, , 1982) 21-22   Source: Authors Thermodynamics: the overall reaction Og We assign arbitrary values to the standard chemical potential for the chemical species variables of the overall reaction Og, taking the following values as reference:   15) is written as Equation ( 16) after using the equation for the chemical potential Here, Q is the reaction quotient of the reaction, and, at the chemical equilibrium, the reaction quotient is equal to the equilibrium constant Q K = , thus obtaining Equation ( 17). ( ) where the standard chemical potential of variable Y is defined from Equation ( 17) as follows: ( ) Finally, for the overall reaction Og, we have the following expressions for the reaction quotient and the Gibbs function, as a function of the reaction extent, i.e., Equations ( 20) and ( 21): Thermodynamics: reaction steps O1 to O5 For the Oregonator reaction steps, we assign arbitrary values to the standard chemical potentials of the variables A and P : The values of the standard chemical potentials for the variables B , X , Y , and Z are assigned according to Equations O1-O5 and ( 17) to guarantee thermodynamic consistency in simulations and calculations: ( ) The overall reaction's thermodynamic equations apply similarly to Equations O1-O5.

Results and discussion
Thermodynamics: reaction steps O1 to O5 Figure 1 shows the dynamics of the variables A , B , Y , and P of the overall reaction Og over time.Curves emerge from the integration of Equation ( 9) to the initial conditions indicated in the legend, subsequently replacing the result in . To better appreciate the result of the integration, we use an equilibrium constant smaller than the Oregonator, with a value of 7 10 .This does not affect the evolution of the Og reaction towards equilibrium.As expected, the species that act as reactants, A , B , Y , decrease until they reach equilibrium, while the product P increases.
From the results of Equation ( 9), and by using Equation ( 19), we compute the thermodynamic quantities describing the evolution of the overall reaction Og towards the equilibrium state.
Figure 2 shows, also based on Equation ( 6), the thermodynamic affinity as a function of the reaction extent for different initial conditions of variable B (i.e., the organic substrate).Note that the thermodynamic affinity decreases continuously and smoothly until it reaches zero at equilibrium, which satisfies the conditions in Equation

Source: Authors
Finally, from the reaction extent and the reaction quotient, we calculate the change in the chemical reaction Gibbs function, r G ∆ , defined in Equation ( 16). Figure 4 shows that r G ∆ increases continuously and smoothly until it becomes equal to zero when the reaction quotient is equal to the equilibrium constant.Thus, from Figures 2, 3, and 4, we verify that the global reaction Og satisfies the requirements of classical thermodynamics in chemical equilibrium, i.e., the thermodynamic affinity as a function of advance is zero, the Gibbs function of the system has a minimum as a function of the reaction quotient, the Gibbs function of the reaction is zero as a function of the reaction quotient, and the equilibrium constant uniquely defines the chemical equilibrium.
The formalism and methodology presented in this section are commonly applied in studying chemical processes.In modeling the homogeneous gas-phase oxidation of sulfur dioxide and titanium (IV) chloride, thermodynamic affinity and the Gibbs function describe the evolution towards chemical equilibrium, approaching zero as a function of the overall reaction extent (Koukkari et al., 2018;Koukkari and Pajarre, 2021), as shown in Figures 2 and 4. A similar approach applies to characterizing the equilibrium of thermal biomass conversion and methanation processes (Kangas, 2015;Kangas et al., 2017).
There are also many applications of thermodynamic formalism devoted to studying the optimization of overall chemical processes in batch reactors, like optimizing the reaction extent and entropy production in sulfuric acid decomposition (Wang et al., 2016).There is also some criticism of the formalism of thermodynamics also.A modified Gibbs function seems necessary to correct inaccuracies in the prediction of equilibrium models, positing that heat exchange between the system and the surroundings is not stoichiometry of the overall reaction Og.This is why the reaction proceeds until the reaction extent is very close to the value of the limit quantity B .

Source: Authors
From the reaction extent obtained using Equation ( 9), we calculate the reaction quotient of the overall reaction Og, defined in Equation ( 20), and the Gibbs function of the system, described in Equation ( 21). Figure 3 shows the

Gibbs function for different initial conditions of variable B
. Note that the Gibbs function decreases as a function of the reaction quotient, which changes until it becomes equal to the value of the equilibrium constant.Thus, the Gibbs function decreases continuously and smoothly until the reaction reaches chemical equilibrium, which is univocally defined when the reaction quotient equals the equilibrium constant.
reversible.As a result, the modified Gibbs function corrects the prediction of the chemical equilibrium and entropy generation of methane steam reforming overall reaction as a function of the reaction extent (Haseli, 2019).A similar approach applies to optimizing ammonia decomposition in hydrogen production (Ojelade and Zaman, 2021).It is important to note that little attention has been paid to including the reaction quotient in the thermodynamic study of chemical processes.The Oregonator: reactions O1 to O5 By integrating the set of Equations ( 10)-( 14), the reaction extent of each step O1-O5 is obtained, and, from these, the change over time of each of the model variables is calculated: This study does not analyze the stability of fixed points in the Oregonator model.Instead, Figure 5 demonstrates that the model displays an oscillatory behavior for the initial variable values provided.Variable B serves as a control parameter for the model's dynamics, decreasing continuously and stepwise until it reaches a constant value, which is consistent with the oscillatory dynamics of the model.The net consumption rate of B depends on which set of reaction steps has a kinetic control on the dynamics: steps O1 to O4 or step O5.The curves in Figure 5 show that, within the constraints established in this study, the model exhibits transient oscillations that eventually disappear, after which the reaction proceeds monotonically.The continuous depletion of B during the reaction process satisfies the second law of thermodynamics, as the primary reactants remain constant, providing the free energy that drives the reaction towards equilibrium.
The variables X , Y , and Z in the Oregonator model exhibit oscillatory dynamics corresponding to intermediate species in the reaction mechanism.X drives the chemical instability in the reaction mechanism, and Figure 6 shows its oscillatory dynamics.The curves (a) and (b) demonstrate that the oscillations in X dampen and eventually vanish, exhibiting monotonic dynamics.These results agree with the well-known experimental evidence of the BZ reaction, which is explained with the FKN mechanism and modeled semi-quantitatively with the Oregonator (Field et al., 1972;Field and Noyes, 1974).Curve (b) highlights the effect of the initial value of B as a control parameter of the reaction system's evolution.The Oregonator transitions from a cycle of damped oscillations with higher amplitude to one with lower amplitude and frequency.This finding is crucial because it indicates that, if the model is not studied for an extended period of time, the completion of an oscillatory cycle could be mistaken for the approach to equilibrium.However, for the constraints imposed in this work, all variables undergo sustained changes, implying that the path to chemical equilibrium is asymptotic and far from coinciding with the end of the oscillatory dynamics.The Oregonator model displays oscillatory dynamics in the variables X , Y , and Z , corresponding to intermediate species in the reaction mechanism.Since each of the Oregonator steps O1-O5 has a reaction extent according to the system of Equations ( 10)-( 14), it is not possible to define the overall extent of the reaction as a function of the model dynamics.To estimate the extent of the Og reaction, we calculated the dynamics of variables A , B , Y , and P , as shown in Figures 5 and 6.The large temporal dynamic of the Oregonator indicates that the extent of the Og reaction calculated on this time scale is also significant and cannot be related to the initial values of the variables.We analyzed the Gibbs function of the Oregonator from the sum of the Gibbs functions of each step O1-O5.Figure 7 shows the variation in the Gibbs function for different initial conditions of variable B as a function of the extent calculated for the overall reaction Og.This result is interesting because it shows how the Gibbs function decreases continuously but not smoothly, describing the oscillatory dynamics of the model.Furthermore, the Gibbs function describes different trajectories in the evolution to equilibrium, depending on the initial conditions.Figure 8 shows how the Gibbs function of the Oregonator reaction, r G ∆ , varies using quantities calculated similarly to that indicated in Figure 7.We obtained r G ∆ by summing the quantities calculated for each step of the model and calculating the reaction quotient from the overall extent of the reaction.The behavior of r G ∆ as a function of the reaction quotient in Figure 8 is like that observed in Figure 4, despite the oscillatory dynamics of the Oregonator.However, even after integrating the model for a long time, r G ∆ does not approach zero, indicating that the system is not at chemical equilibrium.This behavior could be explained by the fact that the estimated reaction quotient has a magnitude that is far from the value of the global equilibrium constant:  8), as shown in Figure 9, where the curves for different initial conditions of variable B demonstrate good agreement between the thermodynamic dissipation of the Oregonator and the reaction quotient.The thermodynamic dissipation decreases continuously and smoothly, approaching zero and indicating equilibrium.However, as in Figure 8, the reaction quotient is not close to the expected value for the equilibrium constant.These results highlight the challenge of applying classical thermodynamics fundamentals to complex reaction systems.While all the thermodynamic fundamentals hold, they cannot be generalized for each step of a reaction system or extended to analyze the overall reaction.As shown in Figures 5 and 6, studying complex processes is challenging from kinetic or energetic approaches, and different methodologies are considered to rationalize the process description (de Oliveira et al., 2016).The traditional method used in this section describes physical, electrochemical, catalytic, and biochemical multistep processes evolving in time (Heimburg, 2021).However, little effort has been made to deepen the thermodynamic study of these systems, since most research focuses on kinetic and dynamic description and analysis (Marin et al., 2019).In the thermodynamics of irreversible processes, entropy production rate is a valuable quantity to optimize the energetic performance of complex processes (Heimburg, 2017;Nieto-Villar, 2020;Arango-Restrepo et al., 2020).Given the formalism used herein, some applications of the entropy production rate show the utility of this quantity to describe multistep chemical systems (Dutt, 1985;Barragán et al., 2015;Barragán and Montoya, 2021).

Conclusions
This study demonstrated that classical thermodynamics principles can be applied to single-step reactions with nonstoichiometric coefficients, as exemplified by the overall reaction Og of the Oregonator model.The thermodynamic affinity decreases continuously and smoothly for a singlestep chemical reaction until reaching zero at equilibrium.In contrast, the Gibbs function decreases until the reaction quotient equals the equilibrium constant for a single-step reaction.
However, for complex coupled reaction systems like the Oregonator, the thermodynamic formalism could not be verified.The lack of knowledge of the reaction extent and the reaction quotient hinders the evaluation of thermodynamic functions and the verification of the state of chemical equilibrium.Non-stoichiometric coefficients in these models, such as parameter f , complicate the determination of the reaction extent because the system's dynamics can remain in non-equilibrium states, distorting thermodynamic quantities.In this context, the rate of entropy production is a valuable quantity.Further research must define and express the reaction extent for coupled chemical reactions.

ϑ=
is the stoichiometric coefficient of the substance i R in the r-th chemical reaction.The i R substances distribute between reactants with a negative ir ϑ coefficient and products with a positive ir ϑ coefficient.The change in the amount of moles i N of the substance i R in the r-th chemical reaction is determined by the reaction extent rThe total change in the amount of the substance i R in the reaction system is calculated via Equation (3):

.
From Equation (4), at constant T and p , the change in the Gibbs function during reaction r is shown in Equation (15). 1 is fundamental in providing consistency to the thermodynamic calculations of chemical equations.For the overall reaction Og in Equation (15), Equation (18) is obtained.

( 7 )Figure 3 .
Figure 3. Gibbs function as a function of the reaction quotient of the overall reaction Og.Curves in the figure correspond to a different initial value for the variable B , as follows: (a)

Figure 1 .Figure 2 .
Figure 1.Dynamics of the overall reaction Og of the Oregonator over time.Curves in the figure are identified with the corresponding variable.The initial values for the variables are 3 0 3 10 , − = × B

Figure 4 .
Figure 4. Change in the Gibbs function of the reaction r G ∆ as a function of the reaction quotient.In the calculations, the following initial conditions are used: 3 0 3 10 , − = × B

Figure 6 ..
Figure 6.Dynamics of the variable X , the autocatalytic chemical species of the Oregonator.Curves correspond to a different initial value for the variable B , as follows: ( ) ( ) 2 2 3.8 10 , 4, 0 10 − − × × a b can calculate the entropy production rate of the Oregonator from the results obtained in Figures7 and 8by using Equation (