Simulation-based investigation of tar formation in after-treatment systems for biomass gasification

Even though biomass gasification remains a promising technology regarding de-centralized sustainable energy supply, its main limitations, namely the issues of unsteady operation, tar formation in after-treatment systems, and consequential high maintenance requirements, have never been fully overcome. In order to tackle the latter two deficiencies and to increase the understanding of thermodynamic and thermokinetic producer gas phase phenomena within the after-treatment zones, a numerical system dynamic model has been created. Thereby, naphthalene has been chosen to represent the behavior of tars. The model has been validated against a wide variety of measured and simulated producer gas compositions. This work particularly focuses on the investigation and minimization of tar formation within after-treatment systems at low pressures and decreasing temperatures. Model-based analysis has led to a range of recommended measures, which could reduce the formation tendency and thus the condensation of tars in those zones. These recommendations are (i) to decrease gas residence time within pipes and producer gas purification devices, (ii) to increase temperatures in low-pressure zones, (iii) to increase hydrogen to carbon ratio, and (iv) to increase oxygen to carbon ratio in the producer gas. Furthermore, the numerical model has been included in the cloud-computing platform KaleidoSim. Thus, a wider range of process parameter combinations could be investigated in reasonable time. Consequentially, a simulation-based sensitivity analysis of producer gas composition with respect to process parameter changes was conducted and the validity basis of the above recommendations was enlarged.


Introduction
The method of biomass gasification via drying, pyrolysis, reduction, consequential partial oxidation to combustible producer gas, gas purification, feeding a gas motor, and producing electricity via a generator has been applied for decades. It remains a promising technology regarding future, decentralized, sustainable energy supply (e.g., Najser et al. [1], French et al. [2], Reed et al. [3], Sansaniwal et al. [4]).
However, the main limitations of wood gasification systems, namely the issue of unsteady operation, excessive tar formation and condensation within after-treatment systems and engines, and consequential high maintenance requirements, have never been fully overcome. In this context, Groenier et al. [5] report that an exemplary 50 kW biomass gasification system requires not less than 30-min maintenance time per day. Reed et al. [3] point out that tar depositions in aftertreatment systems and gas engines are a cause of failure for valves and moving engine parts due to sticking. Asadullah [6] remarks that tar condenses in the low-temperature area in downstream applications, resulting in plugging and fouling of pipes, tubes, and other equipment. Palma [7] highlights that tars frequently cause operational problems to equipment. Furthermore, Palma [7] states that heavy tars (e.g., naphthalene) may condense on cooler surfaces downstream, leading to blockage of particle filters and fuel lines. The author even goes as far as linking the efficient removal of tars from the product gas to the attractiveness of gasification as a whole. And finally, Ruiz et al. [8] declare that the presence of tars in Electronic supplementary material The online version of this article (https://doi.org/10.1007/s13399-020-00915-7) contains supplementary material, which is available to authorized users. the syngas is one of the main technology barriers to the development of gasification.
In order to tackle the challenge of minimizing tar formation as well as to further increase the understanding of thermodynamic and thermokinetic gas phase and tar-related phenomena, this work turns to the development of a system dynamic numerical gas phase model.
Numerical modelling efforts concerning biomass gasification have been undertaken on a wide international basis and for many years. As summarized by Sharma [9], commonly used gasification modelling techniques include the application of thermodynamic equilibrium, chemical kinetics, diffusioncontrolled and diffusion-kinetic approaches, and CFD tools. Boiger [10] extended the range of applied simulation methodology by a system dynamic modelling scheme. The said approach is further elaborated within the current work. With regard to thermodynamic equilibria-related models, Jarungthammachote et al. [11], Ramanan et al. [12], and Roy et al. [13] developed software focusing on downdraft gasifiers. While the former two predicted the composition of syngas in a downdraft waste gasifier for solid waste as well as for cashew shells respectively, the latter highlighted the pyrooxidation zone and kinetically controlled reduction reactions in the reduction zone. Babu et al. [14] developed another equilibrium model in which the effect of oxygen enrichment of air, preheating of air and steam to air ratio on gas composition, reaction temperature, and calorific values were investigated.
By combining the chemical equilibrium and the thermodynamic equilibrium of the global reaction, Melgar et al. [15] predicted the final composition of the producer gas as well as its reaction temperature in a downdraft biomass gasifier. In addition, Sharma [16] used global reduction reactions applying thermodynamic principles based on the stoichiometric approach.
In terms of numerical reaction zone simulation for downdraft gasifiers, Giltrap et al. [17] developed a model for the reduction zone. Babu et al. [18] extended that model by including variable char reactivity. Jayah et al. [19] on the other hand developed a model, which incorporates a flaming pyrolysis sub-model by Milligan et al. [20] along with a gasification zone sub-model. Tinaut et al. [21] described a onedimensional stationary model of biomass gasification in a fixed bed downdraft gasifier. Reed et al. [22] developed a predictive model for stratified downdraft gasification. And Murugan et al. [23] as well as Ngamsidhiphongsa et al. [24] established numerical simulation models for downdraft Imbert gasifiers.
Neither of these models considers tar within the product gas phase. This is either because tar was excluded for the reason of simplification or because it was assumed that in downdraft gasifiers tars are cracked as the producer gas passes through high-temperature zones of the reactor (see, e.g., Sheth et al. [25] and Giltrap et al. [17]). However, Asadullah [6] points out that, even in downdraft gasifiers, tars form during secondary gas phase reactions of devolatilized organic compounds with gasifying agents and generally exist in the producer gas stream. Jaojaruek et al. [26] confirm a lessened but still existing presence of tars within downstream systems of downdraft gasifiers. Authors such as Bhattacharva et al. [27] and Strassen et al. [28] furthermore point to the sensitivity of internal combustion engines towards tar content. While Bhattacharva et al. [27] state that producer gas tar content must be lower than 50 mg/Nm 3 , Strassen et al. [28] and Karuppaswamy et al. [29] mention that up to 100 mg/Nm 3 can be tolerated.
In any case, tars within the producer gas remain an issue to be dealt with, especially if an internal combustion engine is applied further downstream. Thus, any biomass gasification device feeding an internal combustion engine would generally require an appropriate after-treatment system to lower unsatisfactory tar content (see, e.g., Jaojaruek et al. [26]). As a consequence, any simulation model with a focus on producer gas within after-treatment systems should account for the presence of tars as well as for secondary gas phase reactions involving tars. In this context, Barman et al. [30] comment that "as different tars are produced in a gasification reaction through complex set of reactions, predicting tar species in the product gas using any numerical technique is very difficult." As a consequence, the equilibrium-based simulation model presented by Barman et al. [30] involves a representative, rather than a calculated amount of tar, namely 4.5% (mass percentage) of total biomass yield. This value had been chosen on the basis of Yamazaki et al. [31] who investigated the changes in the amount and composition of tar with superficial velocity in a downdraft biomass gasifier. Palma [7] provides an overview of modelling techniques used to investigate tar formation mechanisms within the reactor. As such, the author distinguishes between single compound models, lumped models, and detailed kinetic models.
The model developed within the current work has the goal to better explain, predict, and thus understand naphthalene and tar formation phenomena within after-treatment systems of biomass gasification plants. It advances a system dynamic scheme presented in Boiger [10] by including tars. According to Palma [7], the tar-related aspects within the current model are classified as a "single compound approach." This is because these aspects focus on the tendency of formation and decomposition of one representative tar component: naphthalene.
Being a rather large, aromatic, tertiary tar compound, naphthalene features one of the highest condensation points of all tars to be expected in biomass gasification (see, e.g., Morf et al. [32]). Other authors such as van Paasen et al. [33] also use naphthalene in a similar manner, namely to represent the behavior of un-polar tar species. Additionally, Aldèn et al. [34] found naphthalene to be the predominant condensable tar component even after catalytic cracking. The authors also state that [...] these residual components may cause fouling problems in the downstream equipment of a gasifier.
Morf et al. [32] propose a two-step process to model the formation kinetics of naphthalene out of gravimetric tar via intermediate products such as phenols. In contrast, the herebypresented work omits the demand for accurately modelling the formation and decomposition kinetics of naphthalene. It rather focuses on relative comparisons of thermodynamic tendencies and kinetic rates. Thereby the overall thermodynamic tendency within summed naphthalene formation and decomposition reactions (see, Eqs. [15][16][17][18][19] is interpreted as an indicator for the "likelihood of tar occurrence." Despite these limitations, the chosen approach features the following decisive advantage: while other modelling approaches involving tar can only assume tar content as a fixed value (e.g., Barman et al. [30]), the choice of using naphthalene allows for a dynamic consideration of (representative) tars within the full chemical reacting system. While a previous work by Boiger [35] has focused on modelling conditions within the biomass reactor (similar to Giltrap et al. [17], Jayjah et al. [19], Miligan et al. [20], Murugan et al. [23], and Ngamsidhiphongsa [24]), this work aims at phenomena within the producer gas phase and more specifically within downstream after-treatment systems. Being based on system dynamic principles, the current model includes a combination of calculating (i) divergences in molar Gibbs free energy, (ii) their overall tendency for minimization, (iii) the consequential drive towards chemical equilibrium, and (iv) modelling relative chemical kinetics. The model thus allows for the investigation of both dynamic non-equilibrium and equilibrium states. Thereby, the implementation of kinetic effects relies on Arrhenius laws (as in, e.g., Prakash et al. [36] and Murugan [23]), but ultimately on relative rather than absolute reaction rates (see Boiger [10]).
In order to account for typical process conditions within after-treatment systems, the numerical investigation assumes boundary conditions such as gradually reducing temperatures as well as sudden pressure drops.
Model-based analysis of thermodynamic naphthalene formation tendencies in low-temperature, low-pressure zones has led to a series of results. Namely (i) the statement that system dynamic simulations can qualitatively account for the increase in tar formation at low pressures, which is reported in real-life after-treatment systems (e.g., Palma [7]); (ii) a better general understanding of such effects and ultimately; (iii) the recommendation of a range of constructive and process-based measures to reduce the thermodynamic tendency of occurrence of tars in after-treatment systems.
Furthermore, the latest software version of the system dynamic solver has been adapted to a cloud-based computation environment within the platform KaleidoSim (see Boiger et al. [37]). Thus, massive simultaneous computation of multiple process scenarios has been enabled. On this basis, computation time consumption to simulate a considerably wider range of process parameter settings was strongly reduced. As a consequence, the hereby-presented recommendations could be put into a broader perspective as compared with the status shown already in Boiger et al. [38].

System dynamic modelling concept
Based on system dynamic principles, the newly developed gas phase model is inspired by a unifying perspective on general physical phenomena. The concept thus basically leads back to Gibbs fundamental equation (Eq. 1). Thereby, a change in global energy of any system dE is attributed to the sum of any physical potentials φ j times the change of attributed conservative quantities dΨ j (see, e.g., Job et al. [39]). In the face of thermo-chemical investigation, the global energy change becomes a change of total Gibbs free energy per chemical species ΔG i , the driving potential specializes to the speciesspecific molar Gibbs free energy G i , and the conservative quantity becomes the total amount of moles per chemical species N i (Eq. 1).
In addition, the local temperature T, the molar gas phase composition x i , the local pressure p, and the standard states Θ of the latter two parameters have an impact on the molar Gibbs free energy of each species as seen in Eq. 2.
The calculation of molar Gibbs free energies of formation in standard state ΔG i Θ calls for the knowledge of molar standard entropies and enthalpies of formation, which in turn requires the implementation of the proper thermodynamic coefficients for each species. The latter were extracted from McBride et al. [40]. A balance of species-specific molar Gibbs free energies with respect to stoichiometric constants γ ji for species i and reaction j is the basis for calculating prevailing molar Gibbs free energies of reaction ΔG R,j according to Eq. 3.
Finally, any modelled chemical reaction j is driven by a non-zero molar Gibbs free energy of reaction.

Thermokinetic modelling concept
The rate of chemical reactions I Nj (mol j/s) is modelled as seen in Eq. 4. It is composed of a combination of a reaction-specific relative rate constant f j (mol i 2 /J), the prevailing Gibbs free energy of reaction, and a reaction rate coefficient k j (1/s). The rate coefficient in turn is based on Arrhenius kinetics as provided in Eq. 5. There A j (1/s) is the reaction-specific pre-exponential factor and Ea j (J/Kmol) the reaction-specific activation energy as retrieved from Murugan et al. [23].
The chosen kinetic modelling approach rather relies on relative inter-reaction-comparison than on absolute reaction rates. Thus, a spectrum of varying f j can be investigated, accounting for the reality within a complex heterogeneous thermo-chemical reactor where the speed of individual reactions can vary dramatically in space and time. A related kinetic study on pseudo-equilibria states in wood gas systems caused by a relative variation of f j was published in Boiger [10]. Figure 1 presents a graphical interpretation as well as a very condensed and generalized version of the underlying system dynamic modelling scheme. It combines thermodynamic and thermokinetic aspects. Hereby, the rectangular containers represent conservative quantities Ψ j , the filling heights of these containers are driving potentials φ j and their cross-sectional areas represent system capacities κ j . In addition, the arrows with solid dots are fluxes of the conservative quantities I Ψj , single dots are information processing units (e.g., equations), and the dashed arcs are conveyors of information. Furthermore, N tot and N i are total and species-specific amounts of molecules respectively.

System dynamic modelling scheme
A mathematical interpretation of the graphic solver concept within Fig. 1 amounts to a series of coupled ordinary 1st-order differential equations according to Eqs. 6 and 7, which combine with the coupling relations within Eqs. 2 and 3.
The differential equations are discretized with respect to time. A linear system of equations is assembled and numerically solved by a four-step Runge-Kutta scheme according Hazewinkel et al. [41]. The prototype model was implemented graphically within the system dynamic simulation software Berkeley Madonna. This platform is well suited to obtain an overall visually aided understanding of complex interdependencies within 0D and 1D multiphysics systems. However, Berkeley Madonna has certain limitations regarding overview and usability when it comes to handling larger more complex models. Thus, as the system dynamic solver evolved its codebasis was transferred to a Matlab environment. Finally, in Fig. 1 Graphic interpretation of the essentials of the underlying modelling scheme where i is any chemical species and j any chemical reaction order to establish compatibility with the online cloud platform KaleidoSim, the whole software was again transferred into Ccode.

Chemical species and reactions
The simulation model considers the following interacting chemical species: coal C(s), water H 2 O(g), methane CH 4 (g), carbon-dioxide CO 2 (g), carbon-monoxide CO(g), hydrogen H 2 (g), and naphthalene C 10 H 8 (g). Thereby, naphthalene is used to represent any tar components. The reason for choosing naphthalene as a representative for a wide range of tarmolecules is its relatively high condensation temperature (see also Section 1). If a gasification system is designed such that occurring naphthalene will not condensate, then it is likely that other tar components with lower condensation temperatures will not condensate either.
Those mechanisms occur in combination with (iv) a variety of possible oxidation reactions, which can be summarized as seen in Eq. 11. There, k, n, and m are the relative stoichiometric amounts of carbon, oxygen, and hydrogen atoms respectively.

Solver validation
The solver has been validated in two phases: (i) Validation phase I considering the comparison to another in-house solver based on the method of Lagrangian multipliers (ii) Validation phase II considering the comparison to published experimental and alternative model results

Validation phase I
Within the first validation phase, the system dynamic solver was compared with another in-house solver, which is based on a more "common" approach to find chemical equilibria by minimization of global Gibbs free energy. The said model thus only considers the equilibrium state which is calculated by using the method Lagrangian multipliers as seen in Shabbar et al. [43]. There, the global Gibbs free energy and deviations in atomic species balances are minimized by the introduction and consequential minimization of a Lagrange function L (see Eq. 20).
The Lagrange function consists of two components: one expresses the sum of all Gibbs free energies of formation ΔG i of all molecular species i, while the other stands for the j atomic species balances. Thereby, λ j and b 0 j are the Lagrangian multiplier and the input rate of atomic species j respectively, a i,j is the number of atoms j per molecular species i, and N i is the total number of molecules per molecular species. Figure 2 compares the results of the Lagrange approach to the hereby-presented system dynamic solver in terms of calculated tar-free homogeneous producer gas equilibria compositions at p = 10 5 Pa, an oxygen to hydrogen ratio of R O/H = 1, and a temperature range of 375 K ≤ T ≤ 1350 K. It clearly shows very good correspondence of the results. Discounting gas compounds with molar fractions below 10 −3 , the maximum encountered relative deviation over the entire relevant temperature spectrum amounts to 3.1% of calculated molar fractions. This particular deviation occurred at T = 1025 K, where the system dynamic model yielded x CH4 = 0.02990 while the Lagrangian equilibrium solver yielded x CH4 = 0.02897.
Additional aspects of the validation procedure within validation phase I have already been published in Boiger [10].

Validation phase II
Within the second validation phase, the system dynamic solver was compared with a rather extensive series of published results regarding producer gas composition, based on either experimental data or on alternative models. Experimental data used for validation came from Jayah et al. [44] who measured the producer gas composition of rubber wood at varying moistures and varying air-fuel ratios (mole air per mole biomass, short: A/F); Pratik et al. [45] who gasified wood waste in a downdraft gasifier at varying air-fuel ratios; Simone et al. [46] who investigated the gasification of wood-saw-dust pellets and sunflower seed pressings; Pedroso et al. [47] who presented the measured gas composition of olive wood, peach wood, and pine wood at varying air-fuel ratios; and Dogru et al. [48] who analyzed the gasification output of sewage slug at specific moistures and at varying airfuel ratios.  [49] who modelled the producer gas composition of rubber wood at varying moistures and varying air-fuel ratios; Ptasinski et al. [50] who simulated the gasification of straw and treated wood at varying airfuel ratios; Roy et al. [13] who modelled the gasification of olive wood, peach wood, and pine wood at varying airfuel ratios; and Barman et al. [30] who compared his model to Jayah et al. [44], Dogru et al. [48], Ptasinski et al. [50], Pedroso et al. [47], and Roy et al. [13].
Tables 1, 2, 3, 4, 5, 6, 7, and 8 sum up all essential results in regard to the validation procedure. In order to recreate individual gasification scenarios presented within the validation cases, boundary conditions of system dynamic solver runs were set such that R O/H and R C/H of a wet gas scenario were always in line with the validation data. Since cited literature in general does not clearly state exact temperatures at which gas compositions were measured or calculated, they had to be approximated in most cases. Assuming gas sampling locations towards the entrance of typical gas purification systems, temperatures would range between 700 and 750 K (see, e.g., Martinez et al. [51]).
In the following, the deviations between the system dynamic solver and results from literature are presented twofold: (i) as maximum absolute deviation (Max. Dev.) in mol% and (ii) in terms of the root-mean-square error (RMS) as defined in Eq. 21. Thereby, x i,SD represents calculated results from the system dynamic solver, x i,lit stands for either experimental or model-based results from literature and N is the number of gas compounds.

RMS
The presented comparisons in Tables 1, 2, 3, 4, 5, 6, 7, and 8 have to be evaluated bearing in mind that the system dynamic solver does not include a full reactor model. Instead it mainly focuses on the producer gas phase in interaction with the coal bed as well as on gas-phase reactions in cool down zones of after-treatment systems. Also, consider the fact that R O/H and R C/H values of published experiments and alternative models have been provided as boundary conditions for comparison runs. In addition, temperatures have been chosen such that realistic conditions at probable producer gas analysis locations within after-treatment systems (namely in the entry Table 1 Results of validation phase II for the gasification of rubber wood at varying moistures and air-fuel ratios (A/F). Producer gas compositions as predicted by the system dynamic solver (SD) are compared with experimental results from Jayah et al. [44] and modelled results from Jarungthammachote et al. [49]. Results   area and always between 700 and 750 K) are reflected. Despite these facts, the level of correspondence between final gas compositions predicted by the system dynamic solver and published values covering a multitude of application cases is still quite remarkable. Maximum encountered absolute deviations to comparable experimental runs amount to 5.05% mol (see Table 4, gasification of sewage slug, predicted amount of hydrogen, according Dogru et al. [48]). Maximum encountered absolute deviations to comparable simulations amount to 10.47% mol (see Table 1, gasification of rubber wood, predicted hydrogen content, according Jarungthammachote et al. [49]). The issue that 10.47% mol might appear as a relatively large value can be discarded since a comparison with the corresponding experiments (see Jayah et al. [44]), yields a much smaller maximum deviation of 3.53% mol .

Results and discussion
Based on the, thus validated, thermo-chemical model, a dimensionless timeline of gas being produced, purified, and sucked towards the engine can be simulated. The solver  qualitatively predicts shifts in species concentration as temperature and pressure decrease along the flow path of the producer gas. Naphthalene content is investigated in particular and is used to draw conclusions about tar content in general. On this basis, the following measures are recommended to minimize the tendency for tar occurrence in low-pressure zones: (i) decrease gas residence time, (ii) increase temperatures within the gas after-treatment system; (iii) increase hydrogen to carbon ratio R H/C , and (iv) oxygen to carbon ratio R O/C in the producer gas. While measures (i) and (ii) require modifications to the plant and process design itself (e.g., smaller pipes in combination with stronger suction or recirculation of thermal energy via heat exchangers), measures (iii) and (iv) can be implemented by either removing coal from the reaction zone or by adding either water or process air to the process. Note that these recommendations can only hold if it can be assured that no other essential process parameters (e.g., temperatures, feedstock bulk parameters) are affected by the implemented changes.

Measures to decrease tendency of tar formation and condensation
In the following, simulation results are presented, which lead to the recommended measures formulated above.

Measure (i): decrease gas residence time
The model does not depict chemical reaction kinetics quantitatively but qualitatively and in relative relation between otherwise comparable cases. Thus, simulations are able to show that lower pressure can enhance the tendency towards naphthalene production according to Eqs. 15-19. Figure 3 demonstrates an exemplary simulation run, where lower pressures lead to faster creation of naphthalene as compared with a high-pressure case. In addition, the simulation run behind Fig. 3 also demonstrates that pressure differences do not necessarily have to yield a change in the final equilibrium composition of naphthalene.
In the low-pressure case (Fig. 3, case ii, red), naphthalene is more readily produced along the gas flow path than in the high-pressure case (Fig. 3, case i, green), while the final equilibrium state remains almost the same in terms of naphthalene content. The exemplary result shown in Fig. 3 is not representative of wider process parameter windows. Still, on the basis of being able to simulate individual outcomes like this, it is fair to assume that process condition windows exist in real life, which lead to similar effects in actual gasification systems. The obvious countermeasure to avoid tar condensation within piping systems in the context of this phenomenon is decreased gas residence time in critical low-pressure zones.

Measure (ii): increase temperatures
The recommendation to reduce tar condensation in lowpressure zones by increasing local temperatures is a relatively trivial measure and does not require simulation but is an obvious and effective method to prevent any condensation. As such, the measure is stated here for the sake of completeness. However, increasing temperatures is just a remedy to prevent Table 7 Results of validation phase II for the gasification of peach wood at a specific air-fuel ratio (A/F). Producer gas compositions as predicted by the system dynamic solver (SD) are compared with experimental results from Pedroso et al. [47] and modelled results from Roy et al. [13] as well as from Barman  existing tars from condensation, while measures (i), (iii), and (iv) are meant to prevent or at least reduce tar formation itself.

Measure (iii): increase hydrogen to carbon ratio R H/C
Selected case studies were made to simulate the effect of varying hydrogen to carbon ratios within the fuel input (e.g., cellulose plus additives) on the naphthalene content of the wood gas. Figure 4 shows the comparison of three simulation runs where producer gas with varying R H/C ratios was exposed to temperature and pressure drop. According to this prediction, the thermodynamic tendency of the formation of naphthalene in equilibrium decreases with increasing R H/C . The case with the lowest R H/C yields the highest amounts of naphthalene in the final equilibrium state. Furthermore, the pressure drop causes a clearly discernable spike in naphthalene content (see t * ≥ 4), which points to the correspondence between calculated results and experimentally observed effects (see, e.g., Asadullah [6]). Figure 5 depicts a broader perspective on expected gas phase compositions including naphthalene content, with increasing R H/C . The study reveals that naphthalene content in equilibrium shows a peak at R H/C ≈ 0.49 and decreases continuously with a further increase of the hydrogen to carbon ratio R H/C > 0.49. The said peak shall hereby be referred to as methane-emergence-threshold since it corresponds with the appearance of methane, which is only present in equilibria compositions featuring R H/C > 0.49. Above this threshold, additional hydrogen is apparently used rather for the formation of additional methane, than for naphthalene. Below the methane-emergence-threshold however, additional hydrogen will favor naphthalene formation.
Based on these results, the authors recommend the increase of R H/C in the gasification reactor well beyond 0.49 in order to reduce the thermodynamic tendency for tar formation in lowtemperature, low-pressure zones. This could be achieved, e.g., by either removing coal, which would otherwise continue to participate in gasification reactions (see Eqs. 8-10) from the reaction zone, or by adding controlled amounts of water. Note that this recommendation can only hold if R H/C can be modified without impacting any other relevant process parameters.

Measure (iv): increase oxygen to carbon ratio R O/C
Selected case studies were also made to simulate the effect of varying oxygen to carbon ratios within the feedstock on the naphthalene content of the producer gas. Figure 6 shows the comparison of three simulation runs where gas with varying R O/C ratios is exposed to temperature and pressure drop. According to this prediction, the thermodynamic tendency of naphthalene formation in equilibrium decreases with increasing R O/C . As with R H/C , the case with the lowest R O/C yields the highest amounts of naphthalene in the final equilibrium state.
Again the pressure drop causes a clearly discernable spike in naphthalene content, just as observed in real-life gasification systems. Figure 7 depicts a broader perspective on expected gas phase compositions including naphthalene content, with increasing R O/C . The study shows that naphthalene content in equilibrium continuously decreases with increasing R O/C . At R O/C ≈ 0.55, methane content reaches a minimum, causing a steeper decline in naphthalene content for R O/C > 0.55.
Based on these results, the authors recommend the increase of R O/C in the gasification reactor as a measure to reduce tar content in low-pressure, low-temperature zones. This could be achieved, e.g., by either removing coal from the reaction zone, by adding more process air, or by adding controlled amounts of water. Note that this recommendation can only hold if R O/C can be modified without impacting any other relevant process parameters.

Discussion of recommendations
At first sight, the recommendations to increase R H/C and R O/C factors cannot be backed up by experimental studies that relate fuel input to tar formation (e.g., Walander et al. [52], Kaupp et al. [53], Reed et al. [3]). For example, data out of said literature does not show a clearly discernable tendency for tar reduction in dry output gas with respect to an increasing water supply. However, at closer inspection, it becomes clear that neither these experiments nor any other published experimental works known to the authors are directly comparable to the findings within this work. The reason is that in real-life scenarios it is hard to vary neither R H/C nor R O/C without impacting other essential process parameters. Without an extensive amount of implemented process control, any shift in atomic species ratios will also cause changes in temperature, pressure drop, consistency of the bulk feedstock, and/or charcoal bed as well as exposure of individual particles to reactive oxygen from process air, etc. In contrast, this work isolates the impact on thermodynamic and thermokinetic gas phase reaction tendencies by means of simulation. As such, the above recommendations to reduce the tendency for tar formation and/or Fig. 4 Simulated relative amount of naphthalene x C10H8 (mol) in producer gas versus dimensionless time t*(−), for three cases where R H/C is set to 0.44 (yellow/red), 0.88 (green), and 1.32 (black) respectively, R O/C = 0.706, prescribed pressure (blue) drops from p = 10 5 Pa to p = 10 3 Pa within dimensionless residence time window 3 ≤ t * ≤ 4 and temperature drops from T = 1200 K to T = 500 K within 1.5 ≤ t * ≤ 2. The higher R H/C , the lower the naphthalene content in the equilibrium state. condensation within after-treatment systems can only hold, as long as only the very parameters in question are modified. It is the authors' opinion that, even though this is hard to achieve in actual gasification/aftertreatment systems, these findings constitute a gain in knowledge and understanding in themselves.

Applying KaleidoSim cloud computation to increase perspective on process parameters
As indicated in Boiger et al. [38], the system dynamic model has been converted into a process optimization tool for after-treatment systems of biomass gasification plants. In this context, the computational requirements of the system dynamic solver constituted a certain issue. One single simulation run requires the computation of 202 unknowns. Because of the stiffness of the numerical problem, time-stepping needs to be extremely fine. Thus, one single simulation run would require between 6.67e6 and 1.00e7 time steps. As a consequence, a total amount of up to 202 × 1.00e7 = 2.02e9 unknowns have to be computed in order to reach thermodynamic equilibrium for one set of simulated process parameters. Considering single-precision floating-point numbers according to the IEEE standard, this equals approximately 8.08 Gb of data per simulation run. On this basis, RAM becomes a limiting factor for the solver since at least 10 Gb should be available for each simulation run. Table 9 shows a brief overview of the amount of produced data as well as the number of individual simulation runs required to complete the studies underlying Figs. 2, 3, 4, 5, 6, 7, 8, and 9. By integration into the cloud-based HPC platform KaleidoSim (see Boiger et al. [37]), any issues regarding hardware-limitations could be overcome. In addition, massive simultaneous cloud computing (MSCC) was enabled. Thus, dozens of simulation runs could be executed and evaluated simultaneously and within relatively short time periods. On this basis, a larger perspective on tar minimization measures and their relation to wider process parameter ranges could be gained.  Simulated producer gas equilibrium composition x j (−) with the molar fraction of naphthalene x C10H8 (−) being highlighted on the right axes and with x H2 (−) and x H2O (−) being omitted because they are below 1e−3, versus ratio of hydrogen to carbon R H/C . In all cases, ratio of oxygen to carbon is R O/C = 0.706, prescribed pressure drops from p = 10 5 Pa to p = 10 3 Pa within dimensionless residence time window 3 ≤ t * ≤ 4 and temperature drops from T = 1200 K to T = 500 K within 1.5 ≤ t * ≤ 2 before gas equilibrium is reached. a more general basis for conclusions as well as recommendations for how to run the gasification process with respect to adjusting R H/C and R O/C , the investigated process parameter window was increased. Simulation runs for 0.11 < R O/C < 1.65 were performed where R H/C was parameterized between 0.55 and 1.87. Figure 8 shows the results in terms of predicted naphthalene content within the product gas. This extended investigation confirms that for R H/C > 0.55 the tendency to produce naphthalene decreases for increasing R H/C . This seems to hold true throughout the entire R O/C process parameter spectrum. Thus measure (iii) to reduce tar formation in producer gas systems (see Section 3.1) can be confirmed as stated. On this basis, the statement that increasing R O/C will always lead to a reduced tendency for naphthalene formation (see measure (iv), Section 3.1), must be refined. According to Fig. 8, certain process parameter constellations can be found where a positive gradient of naphthalene content with respect to R O/C appears to exist. Thus, for R H/C ≥ 1.30 and R O/C < 0.8 an increase in R O/C would actually lead to more tars. Outside this relatively narrow process parameter window, higher R O/C would always lead to less tar in the final output gas.

Sensitivity analysis of wood gas composition
In order to compare the relative effect of process parameter variations on producer gas composition in aftertreatment systems, a simulation-based sensitivity analysis has been performed. Condensing the results, a numerically derived, normalized wood gas sensitivity matrix J has been introduced. Basically related to a "standard Jacobian," its entries J ji consist of normalized gas equilibrium compounds x j (P)/x j (P 0 ) being successively derived with respect to process parameters P i . Thereby, x j are the components of the gas composition vector x defined as x CH4 , x CO , x CO2 , x H2O , x H2 , and x C10H8 respectively and P i are the components of the process parameter vector P. The latter components are hereby defined as reactor temperature T R , product gas temperature T P , product gas pressure p, R O/C , and R H/C respectively. On this basis, Eq. 22 shows the full, normalized producer gas sensitivity matrix J. Thereby, P 0 is the process parameter vector at which the sensitivity matrix is to be evaluated. Its components P i,0 are defined as T R,0 = 1200 K, T P,0 = 500 K, p 0 = 10 3 Pa, R O/C,0 = 0.706, R H/C,0 = 0.706 respectively. The gas compounds Fig. 6 Simulated relative amount of naphthalene x C10H8 (−) in producer gas versus dimensionless time t*(−), for three cases where R O/C is set to 0.44 (yellow/red), 0.88 (green), and 1.32 (black) respectively, R H/C = 0.706, prescribed pressure (blue) drops from p = 10 5 Pa to p = 10 3 Pa within dimensionless residence time window 3 ≤ t * ≤ 4 and temperature drops from T = 1200 K to T = 500 K within 1.5 ≤ t * ≤ 2. The higher R O/C is, the lower the naphthalene content in the equilibrium state.
at those process parameters are x j,0 such that x j,0 = x j (P 0 ).
The individual matrix entries J ji in Eq. 22 can be approximated by numerical difference terms ΔJ ji applying a central differencing scheme according to Eq. 23.
ΔJ ji ¼ 1 Fig. 7 Simulated producer gas equilibrium composition x j (−) with the molar fraction of naphthalene x C10H8 (−) being highlighted on the right axes and with x H2 (−) and x H2O (−) being omitted because they are below 1e−3, versus ratio of oxygen to carbon R O/C . In all cases, the ratio of hydrogen to carbon is R H/C = 0.706, prescribed pressure drops from p = 10 5 Pa to p = 10 3 Pa within dimensionless residence time window 3 ≤ t * ≤ 4 and temperature drops from T = 1200 K to T = 500 K within 1.5 ≤ t * ≤ 2 before gas equilibrium is reached.
In Eq. 23, Δp i are process parameter variations which are hereby defined as Δp i = 5 × 10°4 × p i,0 . Executing the central differencing scheme, each entry ΔJ ji requires two separate simulation runs yielding x j (p i,0 + Δp i ) and x j (p i,0 − Δp i ). Figure 9 depicts a graphic evaluation of the results of the producer gas sensitivity analysis in terms of J(P 0 ). Entries relating to x H2O and x H2 are thereby omitted because of their extremely low values in the vicinity of P 0 . The full dataset of the entire sensitivity matrix is compiled in Table 10.
The analysis shows that for the process parameter vector P 0 output gas composition generally features the highest relative thermodynamic sensitivity with respect to changes in R O/C and R H/C . The sensitivity towards changes in temperatures and pressure is comparably low. Furthermore, the evaluation of J(P 0 ) shows that naphthalene content is primarily sensitive to changes in R O/C followed by sensitivity to changes in R H/C . Extrapolating naphthalene behavior to tar formation tendencies, this means that shifting oxygen to carbon ratio is the most effective thermodynamic measure to influence tar formation in biomass gasification systems.

Conclusion and outlook
This work has presented insights into the physical, methodical, and numerical principles behind a system dynamic simulation model for the prediction of naphthalene formation as well as producer gas compositions under  Simulating conditions in the wood gas after-treatment system, the prescribed pressure drops from p = 10 5 Pa to p = 10 3 Pa within dimensionless residence time window 3 ≤ t * ≤ 4. The prescribed temperature drops from T = 1200 K to T = 500 K within 1.5 ≤ t * ≤ 2, before producer gas equilibrium is reached.
varying temperatures, pressures, and atomic species ratios. Besides smaller molecular species, the model includes naphthalene in order to represent the thermodynamic tendency of formation and decomposition of tars with hightemperature condensation points. The model's ability to predict producer gas phase composition was validated against a broad spectrum of published modelled as well as experimental results. Bearing in mind certain limitations, the software was used to derive four concrete measures, which would help to minimize the thermodynamic tendency of tar formation and condensation in low-pressure, low-temperature zones of after-treatment systems. Several exemplary simulation case studies were presented to demonstrate the existence of process parameter windows, which would require those very measures to ensure longer maintenance intervals and smoother, tar-reduced operation of any after-treatment system. Furthermore, an instance of the model was implemented within the cloudb a s e d h i g h -p e r f o r m a n c e c o m p u t a t i o n p l a t f o r m KaleidoSim (see Boiger et al. [37]). This enabled the consideration of wider process parameter windows in reasonable computation time. Thus, the general validity of the proposed naphthalene/tar reduction measures could be refined and the relative effectiveness of process parameter variations could be compared by means of a sensitivity analysis.
Under the name BiogasSim and in combination with KaleidoSim, the hereby-presented system dynamic model will be free to use and accessible via internet browser anywhere in the world by mid-2021.
Code availability Code of system dynamic simulation software is not open source. But source code documentation can be provided upon personal request to the corresponding author. Fig. 9 Visualization of J(‾P 0 ) without entries for x H2O and x H2 . Results of sensitivity analysis of producer gas composition towards process parameter changes. Changes in the composition of gas components with respect to process parameter changes normalized by gas composition at base process parameters‾P 0 Table 10 Full dataset of producer gas sensitivity matrix J according Eq. 22, evaluated at P 0 . Where P 1,0 = T R,0 = 1200 K, P 2,0 = T P,0 = 500 K, P 3,0 = p 0 = 10 3 Pa, P 4,0 = R O/C,0 = 0.706, and P 5,0 = R H/C,0 = 0.706 Software: System dynamic simulation software to be freely accessible via www.woodgassim.ch by end of 2020/mid 2021 within pro-bonoproject "Biogassim" by Kaleidosim Technologies AG. Software is not open source though.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.