THEORETICAL MODEL FOR FURTHER DEVELOPMENT OF INTUMESCENT SUBSTANCES TO REMEDIATE SMOLDERING IN WOOD FIBER INSULATION PANELS

Wood fiber insulation boards, as many other wooden materials, are susceptible to smolder. This type of slow and flameless thermal degradation has three upmost important drawbacks. First, smoldering can develop unseen until damages are noticed; second, it does not need any external heat to keep progressing, thus behaving as a self-sustaining process; third, it may shift into flaming combustion. Although wooden insulation materials are very competitive as insulators, its use is not permitted in several countries beyond mid-rise buildings due to smoldering hazard. As measuring of physical parameters is difficult and expensive at high temperatures, the objective of this investigation was to develop a theoretical model that comprises most relevant physical phenomena in order to serve as a supportive tool for further development of fire-retardant substances. The constructed model presents the novelty that it can simulate the self-sustaining smoldering without needing any external radiation heat, but only the self-heating generated by its own exothermic reactions. The model was built based on a program of experimental testing that included thermo-gravimetric analyses and differential-scanning calorimetry, being able to predict particle degradation at different heating rates and oxygen concentrations with errors of about 7,5 %. The adequacy of the model was also compared at the structural scale against a non-standard cone calorimeter test with terminal switching off heat radiation to emulate self-sustaining smoldering, which was used as model validation showing fits of about 23 % in consideration of mass loss, mass-loss rate and temperature profile. A comprehensive sensitivity analysis comprising 60 distinct parameters permitted to thoroughly assess the influence of each model input parameter, which is being presented as a ranking from the most to the less influencing parameters that prevent or foster self-sustaining smoldering. Several unexpected conclusions raised, positioning species’ densities, capacities and reaction activation energies as the most important parameters. To the best knowledge of the authors, this is the first model that can simulate the self-sustaining smoldering of wooden insulation materials, so it is expected to contribute on further development of fire retardant compounds for wooden products.


INTRODUCTION
Smoldering is the central problem hampering the spread of fibers, chips and other wooden particles as raw insulator material for mid and high-rise buildings in several countries. For instance in Germany, the utilization of wood fiber insulation panels (WIP) is limited up to 7-meter-high buildings due to smolder hazard. Smoldering can indeed be dangerous because, as opposed to a conventional fire scenario, it is a flameless and surprisingly slow thermal degradation process and so it can last unnoticed several hours within walls, facades, roofs and other buildings' partitions until either, irreversible structural damages or critical smoke levels inside compartments are reached. Additional hazard also arises due to potential transitions from smoldering to flaming combustion (Steen-Hansen et al. 2018). The process does not need any external heat source like a flame to keep progressing, which conforms the so-called self-sustaining or self-heating process, in where thermal degradation is solely maintained due to internal exothermic reactions. Apart from this smoldering phenomenon, WIP is an outstanding insulator material with a dry thermal effective conductivity ranging typically from 0,036 W/m•K -0,038 W/m•K (Sonderegger 2011, Brombacher et al. 2012 up to 0,08 W/m•K. In addition, WIP can act as a moisture capacitor in buildings due to the inherent hygroscopic properties of wood, as well as adding mass for improving thermal inertia, among other benefits. Remarkable is the ecological footprint of these products too, as the otherwise widely spread insulators are mostly based on petroleum-based materials (Al-Homoud 2005).
Improving smoldering resistance of wooden products has been a relevant scientific problem for decades and as such, it has been extensively investigated. Multiple experiments and fire-retardant treatments have been developed aiming at stopping smoldering or, at least, improving our understanding of the process. Customary fire retardants comprise boron acid (Wang et al. 2004), silicon oxides (Saka and Ueno 1997), tetraisopropoxytitanium (Miyafuji and Saka 1997), among other substances. However, solely relying on experiments to obtain rational conclusions of this phenomenon is not only difficult and constraining, but also very expensive because it is often not feasible to isolate and thus understand the influence of the multiple physical parameters affecting the smoldering process. Additional complications are experienced due to common uncertainties related to physical measurements at high temperatures (Alifanov 2012). The most important certainties gathered from the experiments being that the complete removal of oxygen stops the smoldering process, and some parameters such as density, porosity, particle size and airflow have a clear influence on the speed of the thermal degradation (Kolb et al. 2017).
Elucidation of physical parameters' influence in smoldering processes was also attempted from a theoretical perspective and several models can be found in the literature, as for instance (Moallemi et al. 1993, Das 1999, de Souza Costa and Sandberg 2004, Guindos et al. 2016, Guindos et al. 2018. Those models are typically based on Computational Fluid Dynamics (CFD) methods such that the smoldering simulation behaves according to the conservation of the mass, momentum, energy and the constituent species of the wooden panel. Simpler models based on simplified equilibrium premises have also been published (Ohlemiller 1991). Several of those models have been successfully validated with experiments -typically standard cone calorimeter experiments -serving us to improve our understanding of the process. However, most of the previous models are constructed and validated under the premise that an external heater (heat source) feeds the smoldering process. For most real scenarios however, there is not such a boundary condition, and the wood fibers are consumed without any external heat in a very slow fashion, i.e. it suffers self-sustaining or self-heating smoldering. Some self-sustaining smoldering models can be found in the literature for substances such as non-aqueous phase liquids (NAPLs) embedded in soil (Switzer et al. 2009, Macphee et al. 2012 or wastewater biosolids (Rashwan et al. 2016) in the context of novel waste processing techniques. For insulation materials however, only models related to self-sustained smoldering of polyurethane foams are published (Walther et al. 2000, Yermán et al. 2017, and, to the best of our knowledge, this phenomenon has not been modelled for wooden insulation boards before. The goal of this investigation was to characterize the smoldering of common low-density wood fiber insulation boards with apparent density of about 170 kg/m 3 -240 kg/m 3 and to create a numerical model capable to simulate the self-sustaining smoldering for those products in order to determine the most influencing parameters. The investigation was organized by first measuring the most relevant degradation parameters at small particles -in where transport phenomena are negligible -, which was termed as 0-D experiments. These measurements comprised the thermal reaction kinetics under different atmospheres and heating rates via thermo-gravimetric analyses (TGA) as well as heat capacities at constant pressure and heat flows with differential-scanning calorimetry (DSC). Secondly, a numerical model based on these measurements as well as the relevant conservation equations was upscaled to simulate structural sized panels. Thirdly, the model was validated via 1-D heat transfer experiment conducted under a typical cone calorimeter test procedure except in that heat radiation was switched off to emulate self-sustained smoldering conditions. Fourth, a parametric study to quantify the influence of the distinct physical parameters was performed.

0-D Experiments as empirical base of the model
In order to build a theoretical model of the smoldering process, the first step consists of defining the relevant species and reactions, i.e. the different solid and gas substances, which play an important role in the thermal degradation process and how these combine with each other. This is customarily achieved in complex materials like wood by setting up a model comprising a set of species (components) and chemical reactions, assuming that the reaction rates can be phenomenologycally represented by Arrhenius-type equations, and attempt to fit the model's mass loss and reaction rates to those measured experimentally via TGA. We also adopted such an approach in this investigation. It is however important to note that, the chemical kinetics involved in pyrolysis of wooden materials comprise hundreds of complicated reactions, whose simulation is not feasible for practical purposes (Yeoh and Kuok 2009). Thus, the crucial definition at this stage from a modelling standpoint is to define up to which simplification level the model still can serve to our practical goals.
Hundreds of distinctly simplified wooden pyrolysis models are available in the literature; only some examples are referred here (Ragland et al. 1991, Di Blasi 1993, Grønli 1996, Bryden et al. 2002, Lautenberger 2007, Anca-Couce et al. 2012. In the process of species modelling, our approach was to simplify chemical kinetics as much as possible. But actually, in our experience at self-sustaining modelling, one cannot make such a decision when looking only at small particles (0-D), but rather one needs to resort from small particles to structural samples (1-D) and vice versa. For instance, as detailed in latter sections, we found out that even though late oxidative lignin reactions may not seem relevant in small particles due to slow degradation rates, they are indeed crucial to self-sustaining the smoldering process at larger specimens.
The reactions, solid species and gas species deemed in our theoretical model are presented in the Figure  1. As shown in the figure, five reactions are considered: (1) the evaporation reaction consisting of the moisture evaporation from wet WIP, yielding to dry WIP and water vapor. This is a highly endothermic reaction peaking its conversion rate about 85 °C; (2) the inert pyrolysis reaction triggers the initial stage of the pyrolysis reaction, which develops without oxygen consumption and converts the dried WIP into early char and pyrolysates (gas products), being in this model considered as a slightly endothermic reaction that peaks about 300 ºC. This reaction may be seen as a global reaction that entails the inert degradation of hemicellulose, cellulose and lignin; (3) the early oxidative pyrolysis is similar to the inert pyrolysis except that it consumes oxygen, is normally produced at slightly higher temperatures, peaking reaction rates are higher, is exothermic and does not generate pyrolysates but products. This can be seen as the first or early part of the oxidative pyrolysis consisting of the hemicellulose and cellulose oxidation; (4) the late oxidative pyrolysis is similar to the previous but happens at higher temperatures, at much slower conversion rates, and involves the lignin oxidation. The competitive simultaneousness of all three pyrolysis reactions typically yields to a 'twin peak' (inert pyrolysis plus early oxidative pyrolysis), followed by long plateau (late oxidative pyrolysis) in the mass loss rates measured in TGA experiments, Figure 2; (5) the char oxidation is a highly exothermic, sharp reaction, peaking about 420 °C and occurring just after conversion of the previous slower reaction. It concerns the degradation of late char using oxygen to produce ash and additional gas products. Accurate modelling of such a sharp reaction, see Figure 2, via Arrhenius expression is unfeasible, not only in terms of curve fitting limitations of exponential functions, but also in the sake of numerical problems related to the convergence of computational models. Note also, that in the model the air is defined as a composition of oxygen and nitrogen. This differentiation allows for varying potential oxygen concentrations in the atmosphere. All the reactions involve the production of one gas species, which implies that the solid phase density always decreases, as shrinking develops.

0-D Experiments
Thermal degradation of small particles presents the advantage that one can assume reactions developing without heat and mass transfer mechanisms, since such small particles react as a 'whole'. Therefore, TGA is the common procedure used to calibrate the phenomenological parameters involved in the distinct Arrhenius parameters used for modelling the reactions. The TGA tests were performed on regular European spruce fibers' insulation boards of density about 213 kg/m 3 produced by the Fraunhofer WKI's technical manufacturing lab, Braunschweig, Germany, from local spruce wood using a dry process. The main production steps included the refining of debarked wood chips in the refiner into wood fibers. By varying the distance between the grinding plates, different fiber size distributions could be achieved. The fibers were then dried to a moisture content of approx. 6 % to 10% in a blow-line. The fibers were weighted and mixed with the PMDI binder in a flying shear mixer. The fibers were then first pre-compressed, then pressed in a hot press at core temperatures of 80 ºC to 100 ° C for 15 minutes to form sheets, and left there for at least 16 hours. Finally, test specimens with dimensions of (100 x 100 x 50) mm³ were cut and conditioned for at least 28 days at 20 °C and 65 % relative humidity, further details of the material manufacturing can be found in Kolb et al. 2017. The tests were performed using the TGA/DSC 3+, Mettler Toledo, Columbus, U.S.A., at the facilities of the Fraunhofer WKI Institute for Wood Research, in Braunschweig, Germany. Each sample weighted about 10 mg; a picture of the raw material used for testing is shown in Figure 2 upon pyrolysis initiation about 220 ºC.
Two series of TGA testing experiments were performed: (1) 9 tests in common air conditions by insufflating 10 mL/min and 40 mL/min of O 2 and N 2 , respectively. Those 9 samples were heated at three different rates of 10 K/min, 20 K/min and 30 K/min from 25 ºC up to 1000 ºC. See the average mass loss and mass loss rate of each three repetitions in the Figure 2. The goal of those tests was to calibrate most Arrhenius parameters considering different heating rates; (2) 9 additional tests were performed at a fixed heating rate of 10 K/min but using three different atmospheres including, 55 % O 2 , 20 % O 2 and 0 % O 2 (100 % N2), see the average values in the Figure 3. The goal of these tests was to check whether the model is capable to simulate the reactions at different oxygen concentrations, as well as calibrating the corresponding oxygen Arrhenius parameters, i.e. n O2 (Equation 4). Apart from the twin-peaks followed by plateau and sharp peak already commented, TGA experiments provided important information related to yields of the distinct solid species. Moisture content of the TGA samples was about 5 % of mass, char yields was typically about 20 %, and ash weighted about 4 % in normal air conditions. From Figure 2 and Figure 3, it is also evident that char oxidation only occurs with oxygen and is less sensitive to oxygen concentration than the oxidative pyrolysis. As presented later, the oxidative pyrolysis is actually the dominating very important for the self-sustaining smoldering. However, in reality this reaction occurs very slowly because it depends on oxygen diffusion.
In addition to the TGA experiments, DSC tests were performed on the same facility and using the same raw material. These tests served to firstly measuring the heat capacity of the wood fibers within a range of 20 ºC up to 120 ºC, being 1710 J/kg•K the average value. Heat capacity showed a linear behavior that can be accurately reproduced with the following Equation 1: Where T r is the reference temperature, set to 295,15 K, and n c,w is 1, and c w,0 is the heat capacity at the reference temperature, in this case 1361,8 J/ kg•K. Secondly, DSC tests were also performed to measure heat flow happening at the different reactions, which served to calibrate the heat of the simulated reactions as further described below. Finally, the effective heat transfer coefficient of the wet panel at ambient temperature via common hot transient bridge procedure in the same laboratory, measured 0,072 W/m•K.

0-D Modelling
The generalized pyrolysis open source finite volume method software Gpyro developed by C. Lautenberger (2014), Reax Engineering Inc., at the University of California Berkeley was used for constructing the numerical model. This is a generalized pyrolysis software used for modelling the thermal degradation of a wide range of materials. The measured heat capacity of Equation 1 was used for the wet insulation material, and the heat capacities for the other solid species (at higher temperatures) were calibrated based on the experiments and using starting values reported in the literature (König 2006). The calibration resulted in heat capacities for dry wood, early char, late char and ash of 1100 J/kg•K, 1100 kg•K, 881 kg•K and 881 kg•K, with power coefficients (n c of Equation 1) of -0,5, -0,5, 0,1 and 0,1, respectively.
The two inert reactions, i.e. evaporation and inert pyrolysis, were modelled using the following Equation 2: Were d ω is the volumetric rate of mass destruction of each condensed species, n is the Arrhenius power factor, A is the Arrhenius pre-exponential factor, E is the activation energy, R is the universal gas constant and T is the current temperature. ̅ ρ t=0 is the average bulk density considering all condensed species at the initial condition, and ̅ ρ is the average density at the current time, being both calculated from the relationship between each solid species' mass fraction, Y i , and bulk density, ρ i (Equation 3): were the indices i and M refer to each condensed species and sum of condensed species, respectively. For the three oxidative reactions however, i.e. early and late oxidative pyrolysis plus char oxidation, the following equation was deemed instead (Equation 4): Equation 4 is similar to Equation 3 except that the volumetric destruction rate d ω is weighted according to oxygen mass fraction, Y O2, and corresponding power factor n O2 , which serves for modelling the importance of oxygen in constraining such oxidative reactions.
Because all five reactions of Figure 1 occur simultaneously, exact analytical solutions of the Arrhenius parameters A, E, n and n O2 of Equation 3 and Equation 4 do not exist or require considerable manipulations. It is therefore necessary to calibrate those parameters based on the mass loss and mass loss rate measured in the TGA experiments via optimization algorithms. In this investigation, we calibrated the values of A, E, n according to the TGA experiments in air conditions at different heating rates, while the n O2 parameter was determined using the TGA experiments at constant heating rate but different oxygen concentrations. The optimization procedure consisted in generating trial seeds of a set of Arrhenius parameters via genetic and shuffled complex algorithms (Lautenberger 2014). Such automated procedure was combined with several manual trials. The results of the optimization procedure for different heating rates are shown in Figure 2, and the corresponding Arrhenius parameters are given in Table 1. The quadratic square error of experimental versus numerical curves (considering both mass loss and mass loss rate) revealed an average fit of 7,5 % after optimization procedure. On the other hand, the simulation of mass loss and mass loss rate at different atmospheres is shown in Figure 3. The optimal calibration of the n O2 parameter is also given in Table 1. While the errors of the mass degradation at different atmospheres were larger than 7,5 % using the same parameters as for the normal air conditions, still the model showed the capability to capture the distinct thermal degradation behavior at different oxygen concentrations. This is evidenced by the lack of char oxidation reaction for the case of 0 % oxygen concentration in the mass loss rate graph shown in Figure 3a, as well as the remaining yields of the mass loss in Figure 3b.
The degree of endothermicity or exothermicity of the reactions was modelled with heat sources' equations. Such equations consist of setting the heat source, Q, in accordance to the rate of each precedent condensed species' destruction, d ω , during each reaction, k, multiplied by corresponding heat release or absorbance of the reaction, ∆H (Equation 5): The heat of the reactions was also calibrated based on experiments by using the DSC-measured heat flows, (Figure 4). The results of such optimization procedure are presented in Table 1.

1-D Theoretical model
Once the 0-D model was built and calibrated, the next step consisted of upscaling the model into 1-D considering the distinct transfer mechanisms via conservation equations.

Solid mass conservation
The mass conservation of the solid phase was rationalized assuming that the average density changes, ∂ ̅ ρ/∂ ̅ t, are given as opposed to the formation rate of the total gas species, fg ω , based on that mass is not created or destructed but just transforms from one phase into another (Equation 6): Where the rate of total gas formation is given as the sum of the vapor formation, ω fv , pyrolysates' formation, ω fpgas , gas products' formation, ω fprod , minus the oxygen consumption, ω dO2 (recall the relationships between species and reactions in Figure 1 Here it is important to consider the stoichiometry of the different species to relate gas formation and destruction with the inert (Equation 2) and oxidative reactions (Equation 4). For instance, the formation of vapor is opposed to wet wood destruction, ω dw (calculated with Equation 2), in consideration to the ratio of dry wood to wet wood densities at the initial time, ρ d0 /ρ w0 (Equation 8): The formation of pyrolysis gas is analogous but given by the inert pyrolysis rate (Equation 4), ω dwi , and the relationship of early char to dry wood densities, ρ ech0 /ρ d0 (Equation 9): The products are generated as the sum of early and late oxidative pyrolysis, ω dewo and ω dwo (both being calculated with Equation 4), and char oxidation, ω dchox (Equation 10): Were ρ ash0 /ρ ch0 indicates the ratio of ash to char densities at initial time, ρ ch0 / ρech0 is the ratio of late to early char densities, and the coefficients 1,41 indicate that 1 gram of dry wood or early char oxidizes with 0,41 grams of oxygen, thus increasing accordingly the yield of gas products. Similarly, the coefficient 3 of Equation 10, indicates that 2 grams of oxygen react with 1 gram of late char increasing the yield of gas products during char oxidation reaction. Those stoichiometry coefficients were estimated according to similar values reported in the literature (Lautenberger 2007). On the other hand, the rate of oxygen destruction is precisely given by the indicated stoichiometry of the three previous oxidative reactions such that (Equation 11):

Gas phase conservation
The conservation of the mass in the gas phase is set up by assuming that the difference of gas concentration considering the average porosity, i.e. ∂ ̅ ρ g · ̅ ψ /∂t, is given by the source of total gas formation, minus the gas convection along the direction of the 1-D model, i.e. ∂m/∂z (Equation 12): Here the average porosity is calculated by weighting the porosities of the distinct species according to volume fractions of each species, X i (Equation 13): Where the volume fractions relate to mass fractions as (Equation 14): and the porosities of each species are calculated as the relationship of the current density of each species in each time, ρ i , and the density that each species would have if a raw material was composed by that sole species and there were not any void at all, ρ s,0,I (Equation 15): In our case we set ρ s,0,i to 4000 kg/m 3 for moist and wet wood, 1900 kg/m 3 for early and late char and 500 kg/m 3 for ash, which yield in Equation 15 to porosities' values similar to those reported in the literature (Grønli 1996).

Solid species conservation
The conservation of solid species is simply given by simple formation and destruction balance for each species, (Equation 16): The destruction of each solid species was already given by Equation 2 and Equation 4 as described in previous sections. On the other hand, the solid species' formation is mostly calculated as the complementary of the gas formation parts. For instance, the formation of dry wood is the complementary part of the vapor formation (note the contraposition to Equation 8), (Equation 17): Similarly, the formation of late char from early char (Equation 18): and the formation of ash from late char (Equation 19): The formation of early char is however, composed by both the formation due to inert pyrolysis and early oxidative pyrolysis (Equation 20):

Gas species conservation
The conservation of gas species is similar to that of solid species, except that the change of mass can also be produced by species convection, ∂(m·Y j )/∂z, and species diffusion, ∂(j j )/∂z (Equation 21): Where the indices j are used to refer to each gas species in contraposition to the subindex i used for solid species. In Equation 21, the gas diffusion of each species is calculated as common Fickian diffusion being the gaseous diffusion coefficients calculated according to Chapman-Enskog theory, see further details in (Lautenberger 2014).

Energy conservation in the solid phase
The conservation of energy in the solid phase is written in terms of enthalpy, h, such that the rate of energy change per volume, ∂( ̅ ρ· ̅ h) /∂t, is given by heat conduction, ∂q/∂z, in-depth radiative heat, ∂q r /∂z,, the exchange of heat with the surrounding gas assuming thermal equilibrium between solid and gas phases, m·c pg ·∂T/∂z, the total heat sources due to thermal reactions as given by Equation 5, i.e. ΣQ k , and the change of sensible enthalpy due to formation and destruction balance of each solid species, Σ( ω fi -ω di )·h i (Equation 22): Where the heat capacity of gases, c pg , was set to 1100 J/ kg•K (estimated from Lautenberger 2007), and the enthalpy of each species, h i , is calculated by simple integration of the corresponding heat capacity of each solid species at a reference temperature up to the current temperature in each time step (Lautenberger 2014), being the average enthalpy of solid species weighted according to mass fractions (Equation 23): The heat flux is calculated by Fourier law (Equation 24): where the average heat transfer coefficient is also weighted according to volume fractions, analogously to Equation 13, setting the temperature variation of each species' coefficient with a power function similar to The radiative heat flux is calculated from average emissivity, ̅ ε, and radiative absorption coefficient, ̅ κ, being both calculated by weighting the corresponding values of each solid species according to volume fractions as defined in Equation 13, see further details of the in-depth radiative heat flux computation in (Lautenberger 2014).

Gas momentum conservation
Finally, the pressure evolution is obtained by defining that the rate of gas density assuming ideal gas law, ∂( ̅ ψ P ̅ M/RT)/∂t, is given due to variations of pressure following Darcy's law, ∂( ̅ Κ ∂P/v∂z)/∂z, sources due to gas formation, ω fg (Equation 7), minus variation due to gravity forces, g·∂( ̅ Kρ g /ν)/∂z, (Equation 26): Where R is the universal gas constant, P is the pressure, and ̅ M is the average molar mass calculated by weighting molar mass of individual gas species with volume fractions as defined in Equation 13. In our case, the molar mass of vapor, pyrolysates, oxygen, nitrogen and products were fixed to 18 g/mol, 44 g/mol, 32 g/ mol, 28 g/mol and 40 g/mol, respectively. The average permeability, ̅ Κ, is calculated also from volume fractions' weighting (analogous to Equation 13). Many distinct values of permeability are reported in the literature, see e.g. (Grønli 1996, Angst andMalo 2010), but it estimation according to model results is a common practice at high temperatures. In our case, we set the permeability of wet wood, dry wood, early char, late char and ash to 1·10 -12 m 2 , 1·10 -11 m 2 , 1·10 -11 m 2 , 1·10 -7 m 2 and 1·10 -6 m 2 , respectively. Note that the increasing values of the effective permeability do not only arise from solid mass density decreasing, but also due to cracks in char and ash. The gas viscosity, ν, was assumed equal to the gas diffusion coefficient calculated according to Chapman-Enskog theory as defined in (Lautenberger 2014). Finally, g is the gravity acceleration, and the gas density is calculated assuming ideal gas law, that is (Equation 27):

Numerical implementation details
As initial conditions, we considered an ambient temperature of 298,15 K, ambient pressure of 101300 Pa and oxygen mass fraction in air of 0,206 based on measurements. The mesh resolution was about 0,25 mm and reference time step was set to 0,04 seconds. Further details of the numerical discretization and solution procedure via tridiagonal matrix algorithm are shown in Lautenberger (2014). Most of the solver tolerances were set to 1·10 -4 . Computation time of the model for small pieces of about 5 cm in conventional computers lasts about 18 minutes.

1-D Validation experiment
In order to validate the 1-D model, we performed a cone calorimeter experiment of a 100×100×50 mm and 213 kg/m 3 WIP sample, made from the same raw material and manufacturer as the small particles used for 0-D measurements, at the facilities of the IBMB of the TU Braunschweig, in Braunschweig, Germany. The cone calorimeter is an extremely precise experiment in where the heat release, gas products, oxygen consumption and mass loss are measured during the entire thermal degradation of a specimen, which is subjected to a constant and controlled radiation of heat at its upper side, see ISO 5660-1:2002ISO 5660-1: (2002. The sample is isolated at both ends, such that heat transfer mostly occurs in a 1-D fashion from the upper to the lower side of the sample, (Figure 5b). However, the experiment we performed was a non-standard test because, after the initial ignition (which occurs just some seconds after experiment starts), we switched off the radiation source of the cone such that the sample was actually burning by itself, thus, showing self-sustaining smoldering as illustrated in Figure 5b. We switched off the radiation because we wanted to simulate the real conditions in where the smoldering process is not fed by any external heat source but rather by its own exothermic reactions. The specimen suffered a complete thermal degradation, as only a marginal mass of about 4% of the initial sample mass, corresponding to the typical ash yields, remained after testing. The self-smoldering of the sample lasted about two hours, whereas the flame, ranging from 2 cm up to 15 cm in height, was only observed during the first 420 s of the test, which was also evidenced by the CO 2 measurements of the exhaust gases. Because the cone temperature is known, the heat radiation over the sample can be calculated without considering the flame. Furthermore, by knowing the height and time period of the flame, an additional re-radiation heat can be calculated (Yeoh and Kwok 2009). This additional heat due to flaming turned to be from 1,1 up to twice the heat released by the cone in our experiment; see a comparison of the heat radiated in the sample with and without considering the initial flame in the Figure 5d. During the test, we also inserted 4 thermocouples type K at 10 mm, 20 mm, 30 mm and 40 mm in depth in order to measure the temperature profile within the sample, ( Figure  5a). Although the test lasted about 2 hours, during the decay phase, as char conversion into ash progressed, the thermocouples slipped through the residual ash yields due to the small density of the remaining material, which made temperature measurements not reliable, see Figure 5c and Figure 5e. Such slipping of thermocouples was only observed after 25 minutes, hence, only the first (reliable) part of the experiment was used to validate the 1-D numerical model.  Table 2 and the reactions' parameters measured at small particles (Table 1). Summary of model referential parameters used for the 1-D simulation of the cone experiment (in addition to those reactions' parameters already presented in the Table 1 and boundary condition with initial flame of Figure 5d).

RESULTS AND DISCUSSION
The comparison of the model vs. experimental results for temperature profile prediction, mass loss and mass loss rate of the 1-D validation experiment is shown in the Figure 6. The model showed an average fit in temperature profile and mass loss of 23 %, which is obviously larger to those errors observed for the 0-D model (7,5 %). Such error increase is mostly attributed due to the high uncertainty in heat capacities, porosity and transfer parameters such as heat transfer coefficient, gas diffusion and permeability among others. An important quantity of this error is allocated within the first part of the heating, up to 100 °C, where water evaporation occurs. Most errors within this region are assumed to happen due to the enormous complexity of water evaporation and re-condensation cycles within the sample, which is not explicitly considered in our model. On the other hand, the comparison of the mass loss rate is complicated without data treatment due to the noisy nature of the experimental curve. In the Figure 6c we show both the real measurement and the smoothed measurement. As it can be seen, in reality the mass loss happens in a continuous but at the same time pronounced sawtooth-like function, most likely due to the extremely mass loss rate sharp peaks produced by the char oxidation reaction, (Figure 2b). This noisy nature of the experimental mass loss curve makes it difficult to compare it with numerical models that rely on exponential (smooth) functions. The model, however, clearly achieved our main goal, which was to predict the self-sustaining response using no external heat except from that radiation used to ignite the sample at the beginning of the test, i.e. by using the boundary conditions of Figure 5d (considering the initial flame only in the first 420 s). This was achieved by using exactly the same parameters calibrated for the 0-D experiments (Table 1). The self-sustaining capability is evidenced by all three curves shown in Figure 6; the temperature profile does not drop down after a long time without external radiation; the mass loss still preserves a pronounced slope up to the end of the measurement; and the mass loss rate tends to a constant value. The previous models of wooden insulation materials always assumed external heat to keep thermal degradation processes, so, even with the aforementioned errors, the model is assumed adequate to discuss self-sustaining smoldering.

Parametric discussion
The main goal of the parametric discussion consists of quantifying what are the most important key parameters that modify the self-sustaining smoldering behavior of WIP. Due to the large number of parameters, it is also important to make such a discussion by varying parameters one in order to rationalize the influences, which mostly is unfeasible to achieve with experiments. It should be clarified that the interpretation of the results herein presented may be up to some extent be dependent of the precise origin of the wood fibers used, as well as the manufacturing process and binding agents. In order to discuss whether parameters make the smoldering 'easier' or 'harder', we selected as referential response the numerical mass loss rate curve of Figure 6c. More specifically, we defined each influence as the average squared sum root difference between the referential curve and the computed one during the 25 minutes that lasted the test, which resulted on a plain scalar percentage. The mass loss curve is deemed as an adequate referential parameter for the sensitivity analysis because, even though in reality it is noisy, the burning rate is a good overall indicator of the self-sustaining smoldering risk, as higher burning rates release larger self-heat and thus are less prone to self-extinguish (Kolb et al. 2017). In this discussion, we are most focused into ranking the relative influence of each parameter in comparison to the remaining ones, rather than quantifying absolute influences. The reason is that, on the development of more effective fire retardant substances, it is of upmost importance to focus chemical and physical modifications on defined targets, and, it is only latter, that one determines up to which extent each particular modification is feasible or not. In addition, because the sensitivity of the distinct parameters is extremely wide -as further presented below -we changed each parameter by ±5 %. This allowed for preserving numerical convergence in the computational fluids dynamics' model, even for the case of the most sensitive parameters, without the need of changing numerical tolerances, mesh, time step, etc.
The summary of the parametrical study is compiled in the Figure 7 and Figure 8. In total, we analyzed the influence of 60 parameters when varied by +5 % and -5 % of the referential values of Table 1 and Table  2. This summed up a total of 120 simulations. The main results of the sensitivity analysis are presented in the following.

Highest effective parameters in reducing mass loss rate
In the Figure 7a, the parameters that are more effective in reducing the mass loss rate when increased by 5 % are ranked from the most to less influencing. Conversely, in the Figure 7b, the parameters more effective in reducing mass loss when reduced by 5 % are sorted. In global terms, the influence of all parameters can be seen as a reciprocal function of the type -1/x, such that three different classes of parameters can be identified.
The first class, which we entitled as highest effective parameters, corresponds to those shown at the upmost left quarter of Figure 7 and feature mass loss decreases larger than 0,5 %. Here the most important parameters having influences larger than 1,5 % in mass loss rate comprise: the increase of dry WIP density, decrease of wet WIP density and increase of inert pyrolysis activation energy. These first results may indeed be surprising. The increase of dry WIP while holding all other parameters fixed, produces an overall increase of wooden mass which needs to be degraded (even when there is less moisture content), which affects in reducing the mass loss rate. On the other hand, the decrease of wet WIP while keeping all other parameters fixed, also generates more proportion of wood to transform, which also diminishes mass loss rate. Another possible explanation to this behavior may be that, for the levels of referential moisture considered, the heating of wood fibers due to water vapor may be very relevant in increasing mass loss rate. Note that due to the extremely high porosity, conduction mechanism is much less effective. The inert pyrolysis activation energy is also very important in the smoldering behavior. This may be surprising because this is a slightly endothermic reaction, however it is only after this reaction, that all other exothermic reactions start, so the 'first resistance' to conversion seems to be crucial.
Under a lesser level of influence we can list the increase of evaporation activation energy, the increase of dry WIP heat capacity, and increase of late char density as well as the reduction of late oxidative pyrolysis heat, all having reductions larger than 1 % but less than 1,5 % on the mass loss rate. The increase of evaporation activation energy severely hampers the first reaction thus affecting the entire chain of reactions and so reducing mass loss rate. The increase of dry WIP heat capacity increases significantly the heat necessary for thermal degradation therefore reducing mass loss rate. The increase of late char density also affects in reducing mass loss, perhaps due to large energy needs in order to trigger late char oxidation reaction. Diminishing the late oxidative pyrolysis heat release also influences highly in reducing the mass loss rate. This seems obvious as this reaction shows a sharp heat release (Figure 2a).
Among the parameters that reduce the mass loss rate between 0,5 % and 1 %, we can list the increase of early char density, the increase of early oxidative pyrolysis activation energy, the increase of the evaporation heat, the reduction of the ash heat capacity and the reduction of the late char oxidation activation energy. Similarly to the effect of late char density, the increase of early char density is effective in reducing mass loss rate, maybe due to increase in energy needs for triggering the subsequent reactions. The increase of early oxidative pyrolysis reaction activation energy could be a reducing parameter because it hampers the reactivity of the entire set of exothermic reactions. The increase of the evaporation heat clearly reduces mass loss rate due to its endothermicity. More difficult interpretation arises from the reduction of ash heat capacity and late char oxidation activation energy. The reduction of ash heat capacity may be influencing in that it yields to faster cooling during the decay phase thus diminishing mass loss rate. On the other hand, the reduction of char oxidation activation energy may be explained because larger reaction rates may increase oxygen consumption, thus diminishing the rates of the remaining oxidative reactions. Summarizing all these highest effective parameters in global terms, one can state that many of the most influencing parameters have something to do with the densities' proportion of the different solid species, as well as selected activation energies, heat of reactions and heat capacities.

Fairly effective parameters in reducing mass loss rate
On the second left quarter of the plots in Figure 7 one can identify parameters of medium influence that decrease mass loss rate from 0,5 % to virtually 0 %. We entitled these as fairly effective parameters. On this region we can highlight the decrease of: heat of the inert pyrolysis reaction, wet WIP heat capacity, ash emissivity, heat capacity of the gas phase, thermal conductivity of the ash, and the way (power) the thermal conductivity of early char changes with temperature. In addition, the parameters that produce medium mass loss reduction when decreased are the heat of the late char oxidation reaction, the Arrhenius power factor for the late oxidative pyrolysis reaction, the early char thermal conductivity, the way (power) early char changes its heat capacity with temperature, the dry WIP thermal conductivity and the oxygen power Arrhenius factor in the late pyrolysis oxidative reaction. From this set of parameters, it is important to remark that the heat of late char oxidation reaction is as important as the heat of the inert pyrolysis reaction, even when the first releases a much larger amount of heat. The reason for this is that it is very important not only to release less heat from oxidative reactions to reduce smoldering, but also that the heat needed for endothermic reactions becomes larger. This additional 'complication' regarding the energy required at the initial steps of thermal degradation seems thus very important to reduce smoldering. Another important aspect here is the influence of oxygen on the late char oxidation. One may expect this influence being larger, however for the oxygen concentrations of typical air, small changes seem not to be of upmost importance. In overall, it is relevant that many of the medium influencing parameters have something to do with heat capacity and thermal conductivity.

Low influencing parameters in reducing mass loss rate
At the right side of the graphs in Figure 7, one can note as almost the half of all the discussed parameters showed an influence much smaller in comparison to the previously discussed ones. We entitled these as low influencing parameters. In brief, among this set of parameters we can highlight the pre-exponential and Arrhenius power factors of most reactions, the porosity (indirectly given by solid density as described in Equation 15) and early emissivities (from wet and dry WIP). It may be surprising that the thermal conductivity of the wet panel has very little influence suggesting that other heat transfers mechanisms such as convection or radiation may be dominant for these materials. Furthermore, it is interesting that the activation energy of the late oxidative pyrolysis has very little influence in comparison to other reactions' activation energies; this may be because this reaction occurs in a very slow fashion in comparison to the others, and so activation energy may become less relevant.

Asymmetry of influencing parameters
One may expect that, if the mass loss rate decreases when increasing a specific parameter, then in case of parameter decreasing, mass loss rate would proportionally increase. However, this is not completely true. In the Figure 8a and Figure 8b, the parameters are ranked according to mass loss rate increase, as opposed to Figure 7. Analogously to the previous discussion, the ranking of the influences may be seen in global terms as a reciprocal function of the type 1/x such that three different groups including highly, fairly and low influencing parameters may be distinguished. However, the specific order or the ranked parameters slightly changes respect to the precedent discussion. This demonstrates the enormous complexity of the problem, as well as certain dependence to a case-specific conditions. Even when parameters are changed one by one, one single change may trigger a chain of modifications in several reactions and transport mechanisms due to the competitive simultaneousness of many involved phenomena. Therefore, the sensitivity analysis should be interpreted as a supportive tool in aiding the development of more effective fire protection substances, rather than an absolute quantification of each single influence.

CONCLUSIONS
An experimental program in small wood fiber particles allowed for calibrating a numerical model capable to capture not only the entire thermal degradation process at different heating rates, but also emulate the main behavioral features of thermal degradation under different oxygen concentrations, as well as heat flows. It was feasible to emulate such widely distinct behavior by defining a unique set of parameters that cast typical Arrhenius (exponential like) equations. Obtaining fits in the order of 3 % -10 % is feasible at this level, provided that a sufficiently detailed group of reactions is elucidated. It is however necessary to account for an intense program of genetic algorithms and/or other optimization tools in order to achieve such fits in a reasonably amount of time.
The main challenge of this investigation was to upscale such small-particle based model, into a model that can simulate real macroscopic wood insulation pieces that are consumed alone, i.e. without any external heat, as they suffer self-sustaining or self-heating smoldering. In order to achieve this, a set of conservation equations, mainly based on the literature, were formulated. With these equations, one can implement the previously mentioned set of reactions and thermal degradation processes within a context of mass and energy conservation, as well as transportation mechanisms. This additional mathematical framework needs assumptions, additional calibrations and a profound literature review, as many of these properties cannot be measured experimentally with sufficient accuracy or its measurement is too expensive.
The adequacy of the structural sized smoldering model was validated against a special test performed in a cone calorimeter. We chose not to perform the typical cone calorimeter test, but rather switching off the cone heat radiation after about 240 s at beginning of the test, letting so the piece to be degraded by self-sustaining smoldering for about 2 hours. Although thermocouples' measurement was not completely reliable due to slipping through low-density ash, the performance of such a test within a cone experiment framework let us to obtain reliable measurements of the smoldering process in a very controlled ambient for about 25 minutes. This data was used valuable to test our smoldering model which was capable to predict the response with a moderate accuracy of about 23 %, but showing clear capability to simulate the self-sustaining smoldering in real conditions without external heat, which was achieved by using exactly the same 'original' set of measured parameters at small particles. Future models should be able to improve the accuracy of our model. Improvements of future modelling efforts should focus on better emulation of the extremely complex evaporation-recondensation cycling of moisture during thermal degradation of the panels. In our experience, this behavior is very difficult to obtain with current modelling premises but seems to be of upmost importance for these highly porous materials. Improvements also should focus on the emulation of the sharp late char oxidation reaction. It seems very restraining to model such sharp behavior with typical exponential-like equations. It would also be helpful to advance in characterization and measuring techniques of physical parameters at high temperatures. The experimental data in this field is scarce, expensive and often case-dependent, such that modelers often experience difficulties in building these types of models, and must accomplish multiple calibrations as well as take considerable assumptions.
A comprehensive sensitivity analysis of 60 modelling parameters, allowed for ranking the most important influences in modifying the mass loss rate at self-sustaining smoldering of typical wood fiber insulation panels. Although the model responses were not completely symmetrical when increasing and decreasing parameter variations, model results were robust and convergence of the model did not fail at any of the 120 simulations. Surprisingly, the solid species proportions of virgin components and yields seem to be the most influencing parameter, as well as activation energies and heat of reactions, and selected heat capacities. We expect that this and other future models can serve to experimentalists in order to develop more efficient fire protecting substances for wood materials in the future.