Detailed Modeling of Kraft Pulping Chemistry. Delignification

This work introduces a phenomena-based model for delignification in the kraft pulping process. The solubilization of lignin is described as a set of chemical reactions representing the entire chemistry of lignin degradation as well as dissolution of the degraded lignin. For modeling, reaction mechanisms and reactions kinetics derived mainly from the literature were used. Each reaction was simulated separately and then combined for the overall degradation. The model was validated with experimental results from pine wood meal pulping under a wide range of reaction parameters. The experimental data presented a good fit with the model. With the aid of the model, the structure and the amount of wood components, in fibers and black liquor, can be determined at any pulping stage. Several engineering parameters can be computed from the detailed chemical composition of liquor and wood or chemical pulp. These include, e.g., kappa number, brightness, yield, active alkali, effective alkali, sulfidity, and higher heating value. ■ INTRODUCTION The kraft pulping process introduced in 1879 is still the dominant method of producing chemical pulps. The purpose of any chemical pulping process is to facilitate the separation of the individual fibers in the feedstock by the dissolution of lignin. In addition to purely chemical events, these processes comprise a number of morphological and physical phenomena. Nowadays, novel approaches in biorefining, which refers to sustainable processing of biomass into a spectrum of biobased products, open new opportunities to improve the current production methods. The traditional pulping industry has already adopted broader biorefinery concepts, for coproduction of biofuels, biogas, and other new products. Additionally, conversion of lignin into such products as chelating agents, antioxidants, pesticides, vanillin, phenols, fertilizers, adhesives, and formaldehyde-free phenolic resin is of high interest in various industries. Computational simulation models are excellent and cheap supporting tools for process development, and several of them have been devised to optimize the pulping process. The first study of kraft pulping kinetics was presented by Vroom. He implemented the concept “H-factor” that combines the pulping temperature and time as a single variable. The idea is that all pulping processes characterized by the same numerical Hfactor value for time and temperature produce pulp of equivalent residual lignin content. Later, Hatton implemented a model that predicts the kappa number and yield, based on the H-factor and effective alkali (EF) charge. This model is applicable to pulping of thin chips of a variety of wood species. The main disadvantages of Hatton’s model is that it does not take into account the effects of varying sulfidity and liquor-to-wood ratio (L:W). Additionally, Hatton’s model does not predict hardwood pulping as well as it does softwood pulping. Vroom’s and Hatton’s models are currently widely used globally in the pulping industry. Kerr implemented a model where the kappa number is predicted as a function of the H-factor. Smith’s kinetic model is the most complex one and divides wood material into five components: high-reactivity lignin, low-reactivity lignin, cellulose, galactoglucomannan, and arabinoxylan. Also, many other models have been created; for example, models by Kleinert, Saltin, Burazin and McDonough, as well as Gustafson and LeMon and Teder divide the delignification process into three phases (initial, bulk, residual) that are described by their own delignification rate equations. Smith and Williams divided lignin into fast and slowly reacting parts in their modeling approach. In their model, the more reactive lignin was thought to possess β-ether bonds while the less reactive lignin was described to contain γ-ether bonds. Nieminen developed a lignin degradation model based on a simplified Purdue model, containing two lignin subcomponents with different reactivities. Bogren presented a general kraft delignification model, based on continuous distribution of reactivity types; however, the model is limited to softwood and its ability to predict industrial-like pulping processes is not satisfactory. Choi’s Received: April 26, 2020 Revised: June 24, 2020 Accepted: June 24, 2020 Published: June 24, 2020 Article pubs.acs.org/IECR © 2020 American Chemical Society 12977 https://dx.doi.org/10.1021/acs.iecr.0c02110 Ind. Eng. Chem. Res. 2020, 59, 12977−12985 This is an open access article published under a Creative Commons Attribution (CC-BY) License, which permits unrestricted use, distribution and reproduction in any medium, provided the author and source are cited. D ow nl oa de d vi a A A LT O U N IV o n O ct ob er 1 2, 2 02 0 at 0 5: 15 :2 5 (U TC ). Se e ht tp s: //p ub s.a cs .o rg /s ha rin gg ui de lin es fo r o pt io ns o n ho w to le gi tim at el y sh ar e pu bl is he d ar tic le s.


■ INTRODUCTION
The kraft pulping process introduced in 1879 1 is still the dominant method of producing chemical pulps. The purpose of any chemical pulping process is to facilitate the separation of the individual fibers in the feedstock by the dissolution of lignin. 2 In addition to purely chemical events, these processes comprise a number of morphological and physical phenomena. 3,4 Nowadays, novel approaches in biorefining, which refers to sustainable processing of biomass into a spectrum of biobased products, 5 open new opportunities to improve the current production methods. The traditional pulping industry has already adopted broader biorefinery concepts, 6 for coproduction of biofuels, biogas, and other new products. Additionally, conversion of lignin into such products as chelating agents, antioxidants, pesticides, vanillin, phenols, fertilizers, adhesives, and formaldehyde-free phenolic resin is of high interest in various industries. 7−9 Computational simulation models are excellent and cheap supporting tools for process development, and several of them have been devised to optimize the pulping process. The first study of kraft pulping kinetics was presented by Vroom. 10 He implemented the concept "H-factor" that combines the pulping temperature and time as a single variable. The idea is that all pulping processes characterized by the same numerical Hfactor value for time and temperature produce pulp of equivalent residual lignin content. Later, Hatton 11 implemented a model that predicts the kappa number and yield, based on the H-factor and effective alkali (EF) charge. This model is applicable to pulping of thin chips of a variety of wood species. The main disadvantages of Hatton's model is that it does not take into account the effects of varying sulfidity and liquor-to-wood ratio (L:W). 12 Additionally, Hatton's model does not predict hardwood pulping as well as it does softwood pulping. 11 Vroom's and Hatton's models are currently widely used globally in the pulping industry. Kerr 4 implemented a model where the kappa number is predicted as a function of the H-factor. Smith's kinetic model 13 is the most complex one and divides wood material into five components: high-reactivity lignin, low-reactivity lignin, cellulose, galactoglucomannan, and arabinoxylan. Also, many other models have been created; for example, models by Kleinert, 14 Saltin, 15 Burazin and McDonough, 16 as well as Gustafson 17 and LeMon and Teder 18 divide the delignification process into three phases (initial, bulk, residual) that are described by their own delignification rate equations. Smith and Williams 13 divided lignin into fast and slowly reacting parts in their modeling approach. In their model, the more reactive lignin was thought to possess β-ether bonds while the less reactive lignin was described to contain γ-ether bonds. Nieminen 19 developed a lignin degradation model based on a simplified Purdue model, containing two lignin subcomponents with different reactivities. Bogren 20 presented a general kraft delignification model, based on continuous distribution of reactivity types; however, the model is limited to softwood and its ability to predict industrial-like pulping processes is not satisfactory. Choi's work 21,22 presented a multiscale modeling system for the delignification process of hardwood chips. His model was able to regulate the kappa number and porosity to desired values. All models mentioned above have shown adequate prediction results; nonetheless many of them are limited by the type of wood (softwood or hardwood). Moreover, the mentioned models do not provide insight into the reaction mechanisms.
In the current work, we aimed to create a model where the real chemistry and other essential phenomena, taking place in kraft pulping, are implemented in the model at the molecular level. This approach means that the various reactions taking place in kraft pulping need to be included in the model in a comprehensive and accurate way, different from the previous kraft pulping models. With the model, it is possible to determine the structure of the wood components in the fiber and black liquor at any stage of the pulping process. The simulator should work as a tool to help in understanding the kraft pulping process in more detail.
■ MODEL Basic Principles. The overall idea of this work was to use the kinetic parameters from lignin model compound studies and apply them for developing a simulation model. Additional hypothesis was done on how to model lignin dissolution.
The main modeling principle was to consider reactions and mass transfer phenomena of true chemical components (water, hydroxide ion, hydrogen sulfide ion, sodium ion, methyl mercaptan, methanol, etc.) and pseudocomponents in the         cases where using true chemical components was not possible. In order to account for the Donnan effect, a two liquid phase approach was used. The wood fiber constituents were assumed to be present in the fiber-bound liquid. The initial amount of water bound to the fiber wall was set according to the fiber saturation point (FSP) determined to be 0.35 kg of water (kg of oven-dried wood) −1 . 23,24 More detailed descriptions of the Donnan effect and mass transfer modeling can be found in previous publications. 25,26 The modeling of changes in the FSP, caused by the degradation and dissolution of the cell wall polymers, is described in Kuitunen's work. 25 Considering the complexity of lignin structure, it was simplified in the model where each lignin pseudocompound represents a group of real components with related structures and chemical properties. The full range of lignin structures used in the model is presented in Table 1 and Figure 1. The macromolecular aspect of lignin was considered in modeling lignin dissolution.
The monomeric lignin pseudounits were assumed to be linked to each other through β-O-4 linkages, which are the dominant type of linkages in lignin. Phenolic and non-phenolic lignin units react at different rates. Thus, phenolic and nonphenolic units were treated separately in the model. The main fragmentation reactions of lignin in alkaline media are the cleavage of β-aryl ether bonds. 1 The amount of β-aryl ether bonds is about 45−50% in softwood lignin and 60% in hardwood lignin. 27 This reaction leads not only to lignin fragmentation, but also to the simultaneous liberation of a new phenolic hydroxyl group. Cleavage of β-ether bonds therefore promotes lignin solubility in two ways: decreasing the molecular weight and increasing the phenolic hydroxyl group content. Lignin dissolution was modeled as a transfer process of molecules from the insoluble fiber wall to the fiber-bound liquid. 28,29 In the current model, lignin−carbohydrate bonds and lignin condensation reactions were not taken in consideration. The wood composition was estimated using the literature on a similar type of wood 1 and from the experimental results 30−32 (Table 2).
In the current work, properties of Scots pine were used for model development and validation. However, it is assumed that the model should be as applicable for softwood as it is for hardwood pulping simulation. The model idea is that wood components and their content entered depending on the wood source and type. Consequently, a difference in hardwood and softwood could be reflected in the model by introduced data on wood composition. The initial composition of chemical liquor varied depending on the experimental setup. The reaction conditions such as temperature, pressure, and consistency were obtained or calculated according to the experiments. Extremely fast mass transfer was assumed; this is the reason why wood meal studies were used for model validation and parameter optimization. The model of mass transfer in wood chips was developed; however, it is not published yet. Computer regression with Kinfit software 33 was first utilized to obtain the unknown reaction rate parameters.
Lignin Fragmentation Reactions. The list of all lignin reactions included in the model is presented in Table 3. All lignin fragmentation reactions were assumed to proceed similarly in the fiber wall and in the liquor, so the same kinetic parameters and reaction routes apply to lignin pseudocompounds in the liquid phase (dissolved) and in the fiber wall. The stoichiometry and kinetics of the lignin fragmentation reactions were based on the lignin model compound studies. The Arrhenius parameters were applied from the literature (frequency factor was calculated based on published activation energies and reaction rates) or regressed against experimental data in the present study. The confidence limit was 93−97% in the regressed values.
Reactions of Phenolic Lignin Units. The reaction of phenolic β-O-4 structures S3 (G) ( Figure 1, reaction a) begins with the reversible formation of an unstable quinone methide (QM). 3 The quinone methide may undergo several reactions depending on the alkaline environment. In the absence of hydrogen sulfide, the degradation of QM proceeds mainly via elimination of formaldehyde, and as a result, the enol ether S5 (αβconjG) is formed ( Figure 1, reaction a). Activation energy for the phenolic lignin degradation in alkaline conditions, with formation of the αβconjG structure, has been reported by Gierer 34 and Kondo. 35 For precise values, kinetic parameters were fitted with Kinfit software and the obtained values (Table 3, reaction 2) were utilized in the model. In the presence of hydrogen sulfide, the formation of QM is followed by a nucleophilic attack of hydrogen sulfide ion which leads to subsequent fragmentation through an intramolecular S N 2 reaction ( Figure 1, reaction a) yielding the thiirane structure S7 (GγβS − ). Gierer, 3, 34 Kondo, 35 and Miksche 36 have reported activation energy values for this reaction. It is known that other competing reactions also take place: in the presence of nucleophiles from lignin or carbohydrates, condensation products may be formed. 37−39 However, this was not included in the model. Alkali-promoted    reaction involving the cleavage of the β-ether bonds by epoxide formation S14 (GαβO − ) and guaiacol S1 (G) liberation may take place, especially at high alkalinity 3 (Table 3, reaction 3).
Reactions of Non-Phenolic Lignin Units. Alkaline cleavage of the β-O-4 linkage in a non-phenolic structure proceeds through the participation of an ionized αor γ-hydroxyl group. 36 As shown in Figure 1b, the reaction leads to the formation of an oxirane S11 (nphGαβO), which then reacts further to S12 (nphGβNu) depending on the nucleophiles present. This is the most important lignin degradation reaction that takes place in the bulk delignification phase. Two studies, conducted by Miksche 36 and Gierer, 34 were taken into account. Parameters were fitted in Kinfit 33 software, in order to obtain the kinetic parameters presented in Table 3 (reaction  6). Nucleophiles, such as HS − and HO − , can attack and open the formed oxirane structure. This reaction promotes efficient delignification by fragmenting the lignin and by generating new free phenolic hydroxyl groups, which play a crucial role in the overall delignification process. 1,40 In cases where nucleophilic sites in lignin or carbohydrates are present, the reaction will result in condensation products. 1 The non-phenolic lignin units are difficult to degrade; however, the alkaline cleavage of β-ether bonds in nonphenolic units can be facilitated by the presence of αor γcarbonyl functionality. 41 The cleavage of β-O-4 bonds in nonphenolic units with α-carbonyl functionality takes place under ambient conditions. 42 The reaction mechanism is presented in Figure 1c. First, structure S2 (αCOnphG) reacts with HO − and forms structure S9 (αCOnphGβγconj) which further reacts with HS − , forming S10 (αCOnphGβγS) as a result of the cleavage of the β-ether bond. The kinetic parameters for S9 formation were published by Miksche 43 and Gierer. 41 Recently we conducted additional experimental work 42 in order to obtain kinetic parameters for the second step and more precise data on the first step of the reaction (Table 3, reactions 4

and 5, respectively).
Demethylation Reactions. Non-phenolic lignin units also undergo demethylation by hydrogen sulfide and methyl mercaptide ions, forming methyl mercaptan and dimethyl sulfide, respectively ( Figure 1, reaction d). 44−48 The kinetic parameters for the demethylation reaction by hydrogen sulfide and methyl mercaptide ions were obtained by fitting kinetic parameters against experimental data presented in the literature 45,47,48 (Table 3,

reactions 7−11).
Modeling of Lignin Dissolution. After modeling the lignin fragmentation reactions, the macromolecular nature of lignin was taken into account in the modeling of lignin dissolution. When describing the conditions for lignin dissolution, it was assumed that the rate of dissolution correlates with (i) the cleavage of β-aryl ether bonds in non-phenolic lignin units and (ii) liberation of new ionized phenolic units. Previously, Tarvo et al. 49 applied a similar lignin dissolution scheme for modeling the chlorine dioxide delignification stage. In that model, the lignin units promoting the dissolution were lignin's acidic reaction products (muconic acid, etc.).
The applied lignin dissolution schemes and rate equations are summarized in Table 3 (reactions 12−21). The dissolution stoichiometry and the kinetic parameters were fitted using Kinfit software using data from wood meal pulping studies. 20,30,32 The low activation energy values (10−30 kJ mol −1 ) indicate that the dissolution was modeled to be fast in comparison with the slower degradation reactions. It should also be noted that the applied schemes provide a mathematical solution to explain the experimental data rather describe truly the complex dissolution phenomena.
Model Validation. The comparison between experimental and modeling results regarding lignin structures could not be done on a molecular scale, because reliable analytical data for detailed characterization of the residual and dissolved lignin were not available. Consequently, the simulation results were converted into standardized variables such as kappa number and lignin yield, in order to make the evaluation. Tarvo et al. 49 presented a detailed explanation of the kappa number and lignin yield calculations in Flowbat software. The lignin content (weight percent) is calculated as presented in eq 1. The kappa number is computed by taking into account only compounds bound to the solid fiber wall. Equation 2 is utilized to convert the content of unsaturated structures in pulp into the kappa number. The compound-specific oxidation equivalent values utilized in the model are listed in Table 4.

■ RESULTS AND DISCUSSION
For the delignification model development, experimental work from the study by Bogren et al. 20,50 was utilized. Since their original experimental data was not reported in the paper, the experimental data points were generated with Matlab software by using their reported model. 20,51 Further in the text, these reproduced experimental data points will be called "exper-  Industrial & Engineering Chemistry Research pubs.acs.org/IECR Article M) HS − ion concentration, the model overestimated the rate of delignification at higher temperatures (≥139°C) ( Figure  4). One of the functions of the HS − ion in kraft pulping is to prevent condensation of lignin fragments by blocking their reactive sites, such as the benzylic carbon. 27 However, these lignin condensation reactions were not included in the delignification model, which could be the reason for the overestimation of the delignification rate at low HS − concentration.
Our model was targeted to be an exceptional tool for detailed modeling; consequently, it is very important to see the influence of the hydroxide ion concentration on the delignification process, in the case of constant hydrogen sulfide ion concentration. Correlation between the model prediction and the experimental results are good as presented in Figures 5 and 6.
In the second set of experiments, the hydroxyl ion concentration was varied from 0.1 to 0.78 M while the concentration of the HS − ion was kept constant at 0.26 M (Figures 3, 5, and 6). In general, the model predicted the time dependence of delignification properly at the applied temperature and hydroxyl ion concentration levels. The only clear exception seemed to the combination of the highest hydroxyl ion concentration (0.75 M) with the highest temperature (168°C ) where the model produced significantly faster delignification than what was experimentally observed ( Figure 5). On the other hand, a similar discrepancy was also seen with the experimental data of Paananen et al., 30,32 who studied the delignification at 1.55 M hydroxyl ion concentration ( Figure  7). In fact, the simulation reproduced the data of Paananen et al. 30 very well, suggesting that the modeling principles were well justified. In all, taking into account the complexity of    The main lignin reactions occur during the first two phases: initial and bulk. In the initial phase, delignification proceeds mainly due to the reaction of phenolic and carbonyl structures. As result, around 20% of lignin degrades in the initial phase. The main part of lignin is removed in the bulk phase. In the residual phase a very small amount of lignin is dissolved. The main reasons for the retarded delignification in the residual phase are reported to be the presence of alkalinestable native lignin structure, the presence of alkaline-stable covalent linkages between lignin and carbohydrates, and the condensation reactions occurring in lignin. 56 Although the model simulated well the initial and bulk delignification stages, it often underestimated the importance of the residual delignification stage. Even though we aimed at modeling the delignification real chemical reactions, we did not include all possible degradation and condensation reactions of lignin and the conditions for its solubilization. While the residual lignin may still contain significant amounts of degradable β-O-4 ether bonds, chemical bonding to the cell wall polysaccharides and stable interunit linkages formed through condensation reactions may prevent lignin from solubilizing. Other features that we did not include in the model include, e.g., aryl group migration reactions and stereochemical and thermodynamic aspects of lignin degradation and dissolution. 52 Because the model is based on real chemical reactions, it is flexible in monitoring the contents of various chemical structures and technical parameters. In Figure 8, simulated kappa numbers are compared with experimental data points from Paananen et al. 53 The fitting was good, and predicted results were close to the experimental data. Other engineering parameters that could be computed from the detailed chemical composition of the reaction system include, e.g., intrinsic viscosity, brightness, and yield of the pulp, as well as active alkali, effective alkali, sulfidity, and the higher heating value of the black liquor.
Certain reactions, as mentioned earlier, are missing from the model, and inclusion of those could improve the model further. Yet another aspect, which could either promote or prevent lignin dissolution, is the presence of lignin−carbohydrate linkages. 53 Part of the work related to carbohydrate chemistry is presented in another article. 55 Lignin−carbohydrate linkages are not taken into consideration in the current work. Besides lignin fragmentation and ionization, lignin dissolution is also dependent on ionic strength, alkalinity, and dissolved wood components in black liquor. Increase in ionic strength during cooking leads to increased delignification rate with increasing pulp yield. Additionally, dissolved wood components promote the delignification and increase in yield. It is known that phenolic lignin units promote dissolution; however, it is also possible that saccharine acid could increase lignin dissolution. 57 The model could be improved by introduction of certain physicochemical aspects of lignin reactivity. For example, lignin reactions in the fiber wall (immobilized     54 has claimed that kraft delignification of wood cannot be described by simple lignin models because in reality the delignification process is more complex. Moreover, the current hypothesis about lignin dissolution could be somewhat inaccurate. It is possible that lignin dissolution should be modeled based on solubility rather than on chemical reactions. Lignin is a polymer, and the movement of large polymer molecules is spatially restricted. That is why chemical degradation is important. Besides the chemical structure, the physical state of a polymer is also important for its solubility properties. Studies 36,52,54 have shown the effect of erythro and threo isomers of β-O-4 structures on lignin degradation and dissolution. The activation energy of the threo isomer was reported to be always higher than that of its erythro counterpart for all phases of delignification. 52 Based on the description above, it would be beneficial to take into account thermodynamics in the delignification process. However, already at this point, simulation results give a reasonable fitting and good correlation with the experimental data.

■ CONCLUSIONS
This paper introduces a detailed model for delignification in kraft pulping, extending earlier work for chemical pulp bleaching and hot water extraction. In this part, lignin fragmentation and dissolution of the fragments were considered. As result of this work, it is possible to fit the model to experimental data very well at different temperatures and chemical compositions of the cooking liquor. However, our hypothesis on lignin dissolution may require additional improvements. At present, the model is able to predict the chemical composition of the fiber wall, as well as the black liquor composition with good correlation. Moreover, such engineering parameters as the kappa number, effective alkali, active alkali, sulfidity, total alkali, and higher heating value can be computed from the detailed chemical composition of the liquor and wood or chemical pulp.