Fischer–Tropsch Synthesis: Computational Sensitivity Modeling for Series of Cobalt Catalysts

: Nearly a century ago, Fischer and Tropsch discovered a means of synthesizing organic compounds ranging from C 1 to C 70 by reacting carbon monoxide and hydrogen on a catalyst. Fischer–Tropsch synthesis (FTS) is now known as a pseudo-polymerization process taking a mixture of CO as H 2 (also known as syngas) to produce a vast array of hydrocarbons, along with various small amounts of oxygenated materials. Despite the decades spent studying this process, it is still considered a black-box reaction with a mechanism that is still under debate. This investigation sought to improve our understanding by taking data from a series of experimental Fischer–Tropsch synthesis runs to build a computational model. The experimental runs were completed in an isothermal continuous stirred-tank reactor, allowing for comparison across a series of completed catalyst tests. Similar catalytic recipes were chosen so that conditional comparisons of pressure, temperature, SV, and CO / H 2 could be made. Further, results from the output of the reactor that included the deviations in product selectivity, especially that of methane and CO 2 , were considered. Cobalt was chosen for these exams for its industrial relevance and respectfully clean process as it does not intrinsically undergo the water–gas shift (WGS). The primary focus of this manuscript was to compare runs using cobalt-based catalysts that varied in two oxide catalyst supports. The results were obtained by creating two di ﬀ erential equations, one for H 2 and one for CO, in terms of products or groups of products. These were analyzed using sensitivity analysis (SA) to determine the products or groups that impact the model the most. The results revealed a signiﬁcant di ﬀ erence in sensitivity between the two catalyst–support combinations. When the model equations for H 2 and CO were split, the results indicated that the CO equation was signiﬁcantly more sensitive to CO 2 production than the H 2 equation.


Introduction
Fischer-Tropsch synthesis (FTS) is a process used to produce a vast range of organic compounds. Discovered nearly a century ago, FTS reacts H 2 and CO over a metal catalyst to produce C 1 -C 70 -chain compounds, including gasoline and alternate materials [1,2]. Thus, it is often used to generate fuels, monomers for common polymers, rubber, etc. Despite the long history of research and industrial use of FTS, it is still unknown how this process works on a microscopic level [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22]. There are at least three proposed mechanisms, namely, the carbide originally proposed by Fischer and Tropsch, CO insertion, and formate mechanisms. Yet, none of these can fully describe the entire FTS process, mainly because of the diversity of materials generated by different active metals [1].
To further comprehend FTS, many have attempted to use mathematical modeling as a useful means of describing the process [23][24][25][26][27][28]. Density functional theory (DFT) calculations, paired with experimental procedures, have provided further insight into the FTS mechanism [18][19][20][21]. These models generally fall into one of two categories-probabilistic models and kinetic models. Probabilistic models attempt to describe the propagation of the carbon chain to generate the products observed [23][24][25]. Kinetic models focus more specifically on the rate of the reaction, determining a mechanism to describe the reaction as a whole [26][27][28]. While these models have significant application for understanding FTS, they do not exist without limitation. Kinetic models depend on the mechanism of the reaction, which, as stated earlier, is not fully understood due to the diverse array of products. Furthermore, they become quite mathematically complex. Probabilistic models focus on the likelihood of a specific carbon chain propagating or terminating. This type of modeling is limited by the monomer used for propagation, which is not fully understood.
This research sought to create a model that describes the rate of change of hydrogen and carbon monoxide as the reaction progresses. Using a series of industrially relevant cobalt-catalyst runs with varying supports, a system of differential equations was created to model the unreacted syngas in terms of the products generated. Each chemical product output term (such as moles of methane, water, carbon dioxide, etc.) was scaled by an impact parameter, a value in units of hr-1. The model was then tested using sensitivity analysis (SA) to determine which impact parameter and, consequently, which associated chemical product term altered the solution of the model the most. From the results of this analysis, conclusions were drawn about the chemical implications, which may then be used to further examine catalytic runs.

Computational Modeling
The computational model was established using chemical engineering process principles from the stirred-tank slurry reactor and basic chemical kinetics. The rate of H 2 and CO coming out of the reactor (unreacted syngas) was modeled using ordinary first-order differential equations (ODEs) [29]. More specifically, the rate of change of H 2 in the reaction can be described by the hydrogen flow into the reactor minus everything coming out of the reactor (H 2 , water, methane, hydrocarbon products). All the products (excluding unreacted hydrogen) were scaled by an "impact parameter". This model was applied to CO as well to produce two ODEs as follows: where Q H2 and Q CO are input rates of H 2 and CO in moles per hour, respectively; q H2 and q CO are the output rates for H 2 and CO in units of hr-1, respectively; MH 2 (t), MCO(t), MH 2 O(t), MCO 2 (t), MCH 4 (t), MC 2 C 4 (t), and Mliq(t) are the time-dependent functions for the moles of H 2 , CO, CO 2 , H 2 O, CH 4 , C 2 -C 4 products, and C 5 + (or liquid) products, respectively; and Z1-10 are the impact parameters in units of hr-1. It is important to note that MH 2 (t) and MCO(t) are the unknown functions for this system. All other mole functions were interpolated from the experimental data (see Supplementary Materials for data), but these two are the functions for which the system is solved. The MH 2 (t) and MCO(t) solutions of the ODE system were compared to the experimentally observed values of these functions (see Supplementary Materials) and the resulting sum of squared error was the "response" for the sensitivity analysis. The * symbol at the end of the hydrogen equation is to note that the bracketed term was added later, thus explaining the variance in sequence of impact parameters. Originally, the carbon dioxide term was not included in the hydrogen equation but was later added to explore the possibility of a water-gas shift (WGS). Due to the innovative nature of the differential equations modeling the reactor, the solver used to evaluate the solutions to this model was the Livermore Solver for Ordinary Differential Equations (LSODA). LSODA performs a test for numerical stiffness as the time steps evolve. A stiff differential equation is one which exhibits solutions with different time scales. For example, stiffness may occur in chemical modeling when simulating reactions with very different reaction rate constants. To solve stiff equations numerically, one must use an implicit solver that requires time-consuming nonlinear equation solutions (e.g., Newton's method for nonlinear systems) for each time step. The LSODA solver primarily uses a non-stiff solver for speed but dynamically checks the solution to determine if a stiff method is required and applies an implicit solver for those time steps. Since the values of the impact parameters were unknown, this approach provided the most flexibility in the numerical solutions. Initial testing of the solution showed good qualitative agreement with the data, indicating that the LSODA solver worked properly with this model. Once the solution was established, sensitivity analysis (SA) was performed. SA is a branch of computational methods that provide insight into how much the total change in the output of a model can be divided up and allocated to the "input" variables/parameters of the model [30]. In other words, SA measures the effect on the output of the model that results from changing the input factors. While there are two types of SA, local and global, local was not used in this work since it involves only measures the change in output at a specified point. Global SA measures the output across a specified range for the input parameters, which was more practical in analyzing this model [30].
While there are many methods for global SA, the Sobol method and Morris method [31][32][33] were of primary interest for application to the model. These are commonly used for such systems and they are both well validated in multiple applications [32,34,35]. The Sobol method is a variance-based global SA method that generates sensitivity indices from ratios of variance terms [33]. The limitation of Sobol is that calculating the variances is very complex mathematically, and the method requires a large number of samples (often more than 10,000) to counteract any error from negative indices. Furthermore, the primary function of the Sobol method is to analyze the interactions by looking at the second-and higher-order indices, but since the primary focus of this work was the impact of a factor on the model output, the Morris method was used as the SA tool for this preliminary study.
The Morris method is known as a one-at-a-time (OAT) global SA method [31]. This method measures the impact of each input factor on the model by changing only one input factor at a time and observing the change in the model output. Specifically, Morris generates a set of elementary effects (EEs) for each input factor across the given range [31]. An EE is generated by starting with the model evaluated at a specific set of input factors. Then, the model is evaluated again with only one factor increased or decreased by a specified amount, D. The unchanged model output is subtracted from the altered model output, the absolute value is taken for the difference, and this absolute value is divided by D, as follows: In the equation above, Y represents the model, X 1 . . . i is the set of parameters, and i = {1, 2, . . . , n}, n being the number of parameters in the model. A set of these EEs using multiple values for D is generated for each input parameter across the given range. Since each set, known as the distribution of EEs, can be quite large, they are sampled. The original Morris sampling technique was a type of random sampling that did not use a high number of samples (10-50), but Campolongo et al. discovered a better sampling technique [31]. In this new technique, 500-1000 original samples are taken from the distribution of EEs. These are then analyzed to determine the optimal 10-50 samples (those that most closely describe the set as a whole). The optimal samples, or optimal trajectories, create a new set that is much smaller and easier to analyze statistically. The average, µ*, and the standard deviation, σ, are taken for each optimal sample set. µ* represents the impact of the specific parameter on the model, and σ represents the interaction between a parameter and the other input factors [31]. The greater the magnitude for these values, the greater the impact (µ*) and interaction (σ) for the corresponding input parameter. Since Morris is a relatively simple method that focuses more on the impact and less on the interactions between parameters, it was chosen as the global SA method for this work. It is important to note that, when performing SA on this model, the impact parameters were the input parameters. While each impact parameter has a corresponding output product (e.g., k 1 corresponds to the hydrogen term for H 2 O output), there is not a necessary relationship between the degree of impact and the product itself. However, this work hypothesizes that there is at least a small direct correlation between the magnitude of impact for an impact parameter and the amount of corresponding product produced.

Results and Discussion
To ensure the effectiveness of this modeling, the mass closure for these runs had to be accounted for. Since this work was to understand, specifically, how the reactants affected each product through the sensitivity analysis, our entire process would collapse and be non-functional if mass closure was not accounted for. To do this, the mass flow controllers were calibrated beforehand, the flow in and out was measured regularly, and all liquids pulled from the system were weighed and accounted for. The difference between output and input displayed almost 96% closure for mass balance for all our work.
Using the data collected from the experimental runs, the mole functions were interpolated, the differential equations were solved, and the solution was then tested using Morris method SA. Initially, the DEs were loosely coupled, meaning that the solutions for MH 2 (t) and MCO(t) could be somewhat dependent on one another and the input parameters (Z1-Z9) could potentially interact. Each set of experimental data was analyzed nine times with Morris, on varying ranges from 0.1-50 to 0.1-500. The number of samples generated varied between 500 and 1000, and the optimal trajectories varied between 0 and 50, as described by Campolongo et al. [31]. The µ* and σ values were collected from the Morris output and the average for all nine analyses was taken. The µ* and σ averages for each data set were compared to the other data sets with the same catalyst-support combination and then compared to the values for the other category of data to discern any variance between the two types of supports.
Next, the two DEs were uncoupled. The CO 2 term for the hydrogen equation, along with k 10 , was added to explore the possibility of WGS occurring in the 0.5 wt% Pt, 25 wt% Co/Al 2 O 3 catalyst. Both equations were analyzed using Morris, following the same procedure as described above. The µ* and σ values were recorded for both the H 2 and the CO equations. These were then compared, primarily the µ* values for k 5 and k 10 , the CO 2 impact parameters for the CO and H 2 , respectively. In theory, should µ* for k 10 consistently be greater than that of k5, it could be indicative of hydrogen having a more substantial role in producing CO 2 . This could only happen through WGS, and thus if Z10 > Z5, it could indicate the possibility of WGS. Otherwise, if Z5 > Z10, this would most likely perpetuate the prevailing theory, primarily for cobalt catalyst FTS, that the Boudouard reaction is producing the CO 2 [36,37].
Listed below are the data collected from the outlined modeling and sensitivity analysis procedure. As a key, Table 1 is the list of the five product categories and the corresponding impact parameters for the two ODEs in the model. For example, Z 4 is the impact parameter for water in the carbon monoxide equation.
Recall that the µ* value for each parameter indicates how much that parameter influences the model. In other words, changing the impact parameter with the highest µ* value will alter the output of the model more than any of the other parameters. The µ* value does not indicate any level of interaction among the parameters, only how significantly that parameter affects the model. To examine the interactions between the parameters effectively, Sobol method SA should be used. Tables 2 and 3 show the data collected from Morris method SA performed on the coupled DEs.
The * is to note that Z 10 was added during the second set of sensitivity analysis runs with the split equations. It does not appear in the first set of data. The last column is an average of µ* across all three runs for each impact parameter. The last column is an average of µ* across all three runs for each impact parameter.
From these data, the order of impact for the nine impact parameters can be observed for both supports. For 20 wt% Co/SiO 2 , the order of impact is Z 7 > Z 3 > Z 4 > Z 1 > Z 6 > Z 2 > Z 8 > Z 9 > Z 5 . In contrast, the order for 0.5 wt% Pt, 25 wt% Co/Al 2 O 3 is Z 7 > Z 3 > Z 4 > Z 1 > Z 6 > Z 2 > Z 8 > Z 5 > Z 9. Although a subtle difference, the order of k 9 and k 5 is inverted, indicating that k 5 had a greater impact on the model with the cobalt-alumina catalyst than the cobalt-silica. These orders of impact held constant across all three sets of data for each type of support, as can be seen in Tables 2 and 3. Table 4 lists the averages for each category side-by-side in order to compare the two supports more closely. From the hypothesis, this variation in the order of the parameter impact indicates that CO 2 , the corresponding product for Z 5 , is more impactful on the alumina support than the silica. Furthermore, based on the proposed correlation between the products and corresponding parameters, alumina produces a greater amount of carbon dioxide than does silica. To verify this hypothesis, the original data from the data sets were examined, and the theory was confirmed. As Table 5 shows, the cobalt-alumina catalysts did in fact produce consistently more carbon dioxide than the cobalt-silica. Due to this CO 2 production difference, the two DEs were separated and examined individually using the Morris method. Only the Category 2 data were explored since the CO 2 parameter was more impactful for the cobalt-alumina support. The additional carbon dioxide must come from an additional source. According to the hypothesis, that source could be either water-gas shift (WGS) or the Boudouard reaction. In order to account for any CO 2 produced as a result of hydrogen interaction, thus allowing for the possibility of WGS, the parameter Z 10 was added to the hydrogen equation along with the CO 2 mole function. Tables 6 and 7 show the data collected from Morris method SA for the separated equations. Hydrogen equation only. The last column is an average of µ* across all three runs for each parameter. The row in bold is the parameter of interest (i.e., the parameter for CO 2 ).
The parameters Z 5 and Z 10 were of primary interest. As is evident from the data, the impact of k 5 was invariably greater in magnitude for the CO equation than the impact of Z 10 for the H 2 equation across all three data sets in Category 2. Since Z 5 > Z 10 , the likelihood of carbon dioxide being produced from WGS is very minimal. The results indicate that the extra CO 2 is a product of the Boudouard reaction. Carbon monoxide equation only. The last column is an average of µ* across all three runs for each parameter. The row in bold is the parameter of interest (i.e., the parameter for CO 2 ).

Catalyst Preparation
The cobalt and platinum starting materials were cobalt and platinum nitrate purchased from Alfa Aesar (Haverhill, MA, USA). A 25% cobalt alumina catalyst was prepared through a slurry impregnation procedure using Catalox 150 (high purity γ-alumina, ≈150 m 2 /g, Sasol, Hamburg, Germany) as a support. Next, the calcination of the catalysts was performed under air at 623 K for 4 h using a tube furnace (Lindberg/MPH, Riverside Road Riverside, MI 49084). The alumina was calcined for 10 h before impregnation and then cooled under an inert gas to room temperature. Co(NO 3 ) 2 ·6H 2 O (99.9% purity) was used as the precursor for Co. In this method, which follows a Sasol patent [38], the ratio of the volume of solution used to the weight of alumina was 1:1, such that approximately 2.5 times the pore volume of solution was used to prepare the loading solution. A two-step impregnation method was used allowing a loading of two portions of 12.5% Co by weight. After the impregnation, the catalyst was then dried under vacuum by means of a rotary evaporator (Buchi R-100, Flawil, Switzerland) at 353 K, then slowly increased to 368 K. After the second impregnation/drying step, the catalyst was calcined under air flow at 623 K for 4 h.
The 20% Co/SiO 2 catalyst was also prepared using the aqueous slurry phase impregnation method, with cobalt nitrate as the cobalt precursor. The support was PQ-SiO 2 CS-2133, surface area about 352 m 2 /g. The catalyst was calcined in flowing air or flowing~5% nitric oxide [39] in nitrogen at a rate of 1 L/min for 4 h at 623 K.

BET Surface Area and Porosity Measurements
Surface area, pore volume, and average pore radius of the catalyst calcined at 623 K was evaluated through Brunauer-Emmett-Teller (BET) by using a Micromeritics Tri-Star 3000 gas adsorption analyzer system (Norcross, GA, USA) ( Table 8). About 0.35 g of the catalyst was first weighed out, then added to a 3/8 in. o.d. sample tube. The adsorption gas was N 2 ; sample analysis was performed at the boiling temperature of liquid nitrogen. Prior to the measurement, the chamber temperature was gradually increased to 433 K, then evacuated for several hours to approximately 6.7 Pa.

Hydrogen Chemisorption
Hydrogen chemisorption was conducted using temperature-programmed desorption (TPD), Table 8, with a Zeton-Altamira AMI-200 instrument (Pittsburgh, PA, USA). The Co/Al 2 O 3 (≈0.22 g) was activated in a flow of 10 cm 3 /min of H 2 mixed with 20 cm 3 /min of argon at 553 K for 10 h and then cooled under flowing H 2 to 373 K. The sample was held at 373 K under flowing argon to remove and/or prevent adsorption of weakly bound species prior to increasing the temperature slowly to 623 K, the temperature at which oxidation of the catalyst was carried out.
The TPD spectrum was integrated and the number of moles of desorbed hydrogen determined by comparing its area to the areas of calibrated hydrogen pulses. The loop volume was first determined by establishing a calibration curve with syringe injections of hydrogen into a helium flow. Dispersion calculations were based on the assumption of a 1:1 H:CO stoichiometric ratio and a spherical cobalt cluster morphology. After TPD of hydrogen, the sample was reoxidized at 623 K using pulses of oxygen. The percentage of reduction was calculated by assuming that the metal reoxidized to Co 3 O 4 . Results of chemisorption are reported in Table 9. Table 9. H 2 chemisorption (temperature-programmed desorption (TPD)) and pulse reoxidation of metallic phases of cobalt-supported catalysts (cat.).

Catalyst Activity Testing
Reaction experiments were conducted using a 1 L stirred-tank slurry reactor (STSR; Fort Worth, TX, USA) equipped with a magnetically driven stirrer with turbine impeller, a gas inlet line, and a vapor outlet line with a stainless steel (SS) fritted filter (2 mm) placed external to the reactor. A tube fitted with a SS fritted filter (0.5 mm opening) extending below the liquid level of the reactor was used to withdraw reactor wax (i.e., wax that is solid at room temperature), thereby maintaining a relatively constant liquid level in the reactor. Separate Brooks Instrument mass flow controllers (Hatfield, PA, USA) were used to control the flow rates of hydrogen and carbon monoxide. Carbon monoxide, prior to use, was passed through a vessel containing lead oxide on alumina to remove traces of iron carbonyls. The gases were premixed in an equalization vessel and fed to the STSR below the stirrer, which was operated at 750 rpm. The reactor temperature was maintained constant (≈274 K) using a temperature controller.
A 12 g amount of the oxide cobalt catalyst (sieved 63-106 µm) was loaded into a fixed-bed reactor for 12 h of ex situ reduction at 623 K and atmospheric pressure using a gas mixture of H 2 /He (60 L/h) with a molar ratio of 1:3. The reduced catalyst was transferred to an already capped 1 L STSR containing 310 g of melted Polywax 3000 (Baker Hughes, Houston, TX) under the protection of inert nitrogen gas. The reactor used for the reduction was weighed three separate times, before and after reduction, and after catalyst transfer, in order to obtain an accurate amount of reduced catalyst that was added. The transferred catalyst was again reduced in situ at 503 K at atmospheric pressure using pure hydrogen (20 L/h) for another 24 h before starting the FT reaction. Each run was held at 493 K and 18.7 atm (275 psig). The CO:H 2 ratio was 1:2 for each experiment. The weight hourly space velocity (WHSV) varied from 3.0 to 6.0 slph/g catalyst, depending on the purpose of the experiment.
Gas, water, oil, light wax, and heavy wax samples were collected daily and analyzed. Heavy wax samples were collected in a 473 K hot trap connected to the filter-containing dip tube. The vapor phase in the region above the reactor slurry passed continuously to the warm (373 K) and then the cold (273 K) traps located external to the reactor. The light wax and water mixture were collected daily from the warm trap and an oil plus water sample from the cold trap. Tail gas from the cold trap was analyzed with an online HP Quad Series Micro gas chromatograph (GC) (Palo Alto, CA, USA), providing molar compositions of C 1 -C 5 olefins and paraffins, as well as H 2 , CO and CO 2 . The analysis of the aqueous phase used an SRI 8610C GC (20720 Earl St. Torrance, CA 90503, USA) with a thermal conductivity detector (TCD).
Products were analyzed with an Agilent 6890 GC (Santa Clara, CA, USA) with a flame ionization detector (FID) and a 60 m DB-5 column. Hydrogen and carbon monoxide conversions were calculated based on GC analysis of the gas products, the gas feed rate, and the gas flow that was measured at the outlet of the reactor.

Conclusions
In conclusion, the mathematical model created in this work accurately described the FTS process in terms of rate of change of syngas. When the model was analyzed using Morris method SA, a significant difference was observed between cobalt-silica and cobalt-platinum-alumina catalysts, namely a difference in the effect of the impact parameters k 5 and k 9 . Based on the initial hypothesis, this difference indicated that the cobalt-platinum-alumina generated more carbon dioxide than did the cobalt-silica. This was verified by the raw data as well. These conclusions are highly specific to the two catalyst-support types in question, and thus more general conclusions must come from further research. Thus, the equations were separated, and the cobalt-platinum-alumina catalyst data were reexamined, now with an additional CO 2 term on the hydrogen equation to test for the possibility of WGS. The observed results showed that the carbon monoxide contribution to CO 2 was consistently higher than the hydrogen contribution. Since WGS is the only pathway in which hydrogen produces CO 2 during FTS, WGS was excluded as a source of CO 2 production. Hence, the Boudouard reaction was determined to be the more probable source of additional carbon dioxide. Though the Boudouard reaction yields a high amount of carbon on the surface, to the best of our knowledge, there is no direct evidence that this reaction leads directly and solely to inert graphitic carbon (coke). Moreover, one of the main mechanistic routes proposed for FT is the carbide mechanism, where CO completely dissociates before reacting with hydrogen. In order to understand this rate of formation for C, evidence of more than one type of carbon on the catalyst surface, an active phase and an unreactive phase, has been observed on multiple catalysts such as cobalt [40] and Fe [41].
Further indications using sensitivity analysis show that both the H 2 and CO equations displaying Z 3 and Z 7 obtained the highest value, indicating the liquids are most sensitive to any deviations brought on by H 2 and CO. Though expected, these details could further indicate the importance for a specific balance in the FTS system between the two reactants, suggesting the difficulty in tuning a specific catalyst for C 5 + products as a whole. Additionally, the C 5 + products are more sensitive in the CO reaction than the H 2 , which implies that the balance for the FTS system lies more heavily in the CO than for H 2 , primarily describing the importance of CO, possibly due to the metal's ability to back donation [42] in the FTS process. While this research is preliminary, the model and analysis produced some excellent descriptive results for FTS. Although this work only examined cobalt catalyst data, future work could include the incorporation of multiple catalysts and supports. Impact parameters could potentially be optimized in future work to find the ideal values for a given catalyst. The model can be adapted to include terms for olefins and paraffins, and the liquid (C 5+ ) term could be broken down into individual terms for each number of carbons per chain (e.g., C 5 , C 6 , C 7 , etc.). The equations could be changed to predict specific product amounts based on initial conditions such as catalyst, support, temperature, pressure, flowrate, etc. Ultimately, this type of FTS modeling not only describes the process in various and versatile ways but can also be adapted and altered to potentially produce a predictive model. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.