Elucidating the characteristic energy balance evolution in applied smouldering systems

Applied smouldering systems are emerging to solve a range of environmental challenges, such as remediation, sludge treatment, off-grid sanitation, and resource recovery. In many cases, these systems use smouldering to drive an efficient waste-to-energy process. While engineers and researchers are making strides in developing these systems, the characteristic energy balance trends have not yet been well-defined. This study addresses this topic and presents a detailed framework to uncover the characteristic energy balance evolution in applied smouldering systems. This work provides new experimental results; a new, validated analytical description of the cooling zone temperature profile at steady-state conditions; insight into the characteristic temperature changes over time; a re-analysis of published data; and a robust framework to contextualize the global energy balance results from applied smouldering systems. Altogether, this study is aimed to support researchers and engineers to better understand smouldering system performance to further the development of environmentally beneficial applications.

Smouldering is a flameless form of combustion driven as oxygen in the pore-space transports to the surface of the fuel and then reacts [46,47].Smouldering can be self-sustaining so long as the energy released rate overwhelms local heat losses rate (e.g., due to heating nearby fuel, phase change processes, endothermic reactions) [46,[48][49][50].Beyond the pore-scale, smouldering is characterized by a series of zones that evolve throughout space and time.Fig. 1 presents a conceptual model of these zones in an applied smouldering system, along with a common centreline temperature profile and a picture of an oil-drum sized reactor (DRUM).These systems often use wastes mixed within inert porous media (e.g., silica sand) to provide the necessary mixture conditions for smouldering, e.g., (i) permeable pathway for oxygen, (ii) sufficiently high specific surface area, and (iii) energy recycling between the sand, air, and fuel [23,46,51].
The zones in Fig. 1 are common to smouldering systems, i.e., the cooling, reaction, and inert heating zones.The reaction zone is the driving force that delineates all chemical reactions, i.e., exothermic reactions (primarily fuel oxidation) and endothermic reactions (primarily pyrolysis) [46].Though these reactions may compete throughout space, pyrolysis generally precedes oxidation [52,53].Ahead of the reactions, the energy evolved from smouldering drives fuel preheating and phase change processes (e.g., boiling/evaporation [54]) in the inert heating zone.Behind the reactions, the remaining hot inert material is cooled by the incoming applied air and perimeter heat losses [29,55].The heat transfer in this zone is fundamentally driven by local thermal non-equilibrium (LTNE) [29,[56][57][58].
While Fig. 1 illustrates there are multiple dynamic processes that evolve in smouldering systems throughout space, it is also useful to consider the energy balance around the entire system, i.e., the global energy balance.Equation (1) presents the key global energy balance terms that govern the waste-to-energy potential of smouldering systems [48,59]: where Ėin and Ėout are the rates of energy added for ignition and lost primarily via convection upon smouldering completion, respectively; Ėoxid and Ėpyr are the rates that energy is chemically released (from exothermic reactions, primarily oxidation) and absorbed (from endothermic reactions, primarily pyrolysis), respectively; Ėloss is the rate of energy lost as perimeter heat losses; and dE net /dt is the net rate of energy accumulation in the system.Note that Ėin and Ėout only represent initialand end-effects, respectively, and can be neglected throughout most of smouldering propagation [27,29].In addition, the pyrolysis effects are often small in robust applied smouldering systems, e.g., using bitumen [48,59].Alternatively, the net effect of Ėpyr and Ėloss can often be  combined into an 'effective' source term in applied smouldering systems where these reactions travel close together [46,60].Therefore, the dE net / dt in applied smouldering systems is often governed by the balance between Ėoxid (or an effective source term) and Ėloss [27,29].Note that Eq. ( 1) does not account for energy lost due to water evaporation, but it could be included [54].
In addition, it is useful to examine the terms from Eq. (1) integrated over time, i.e., E = ∫ t 0 Ėdt.For example, the energy accumulated in smouldering systems, E net , is solved by either balancing all terms in Eq. ( 1) and integrating over time, or integrating temperatures throughout time and space [27,29]: Torero et al. [46], recently provided an overview of smouldering energy balance developments, which is summarized below.The smouldering front characteristics ultimately result from a local energy balance, i.e., only within the reaction zone, which was resolved by many researchers using traditional combustion analyses, e.g., by Dosanjh et al. [61] and Aldushin et al. [62].However, the local energy balance in smouldering systems often relies on dynamic inputs that need to be resolved across the entire system, i.e., using a global energy balance [46].Zanoni et al. [48,59], provided a numerical modelling approach to resolve the terms from Eq. ( 1) both globally and locally.A major advantage of this approach is that it provided insight into the dynamic evolution of smouldering systems and the impacts of system input changes on smouldering performance.For example, Torero et al. [46], provides examples that demonstrate how global and local energy balance results can together predict system quenching and the influence of interventions to maintain self-sustaining smouldering.Multiple recent studies have since extended the insights from Zanoni et al. [48,59], to better understand many phenomena, e.g., (i) limits of ignition/extinction [50,63,64] (ii) sensitivities to key operating parameters (e. g., fuel concentration, and air flux) [10,65], and (iii) system energy efficiency gains with scale [27,29].However, many of these studies have focused on smouldering characteristics that were far from global steady-state conditions, and few studies have pursued steady-state analyses.Fig. 2 presents numerical modelling results from Zanoni et al. [48], which resolved all terms from Eq. (1), as an example to demonstrate key global energy balance changes over time common to applied smouldering systems.Moreover, this figure provides modelling data that simulates smouldering system behaviour in a large domain that facilitated effective global steady-state conditions (note that these model results are re-analyzed in Section 3 as LAB B1 -see Table 1).Fig. 2a shows the temperature-time profiles from the simulation, Fig. 2b shows the evolutions of the terms on the righthand side of Eq. ( 1) resolved over the entire system, and Fig. 2c shows the evolutions of the balance of all terms in Fig. 2b (i.e., rate of net stored energy, black curve) and the net stored energy (i.e., red curve).As mentioned above, Fig. 2b shows that Ėin and Ėout were indeed initial-and end-effects, respectively, and Ėpyr was minimal.Therefore, Ėoxid and Ėloss governed the energy balance.That is, the system initially accumulated large amounts of energy in the 'early, transient times' as Ėoxid > Ėloss .However, while Ėoxid was nearly constant, Ėloss increased over time and drew more energy from the system as the cooling zone grew and provided a longer region for heat losses.Therefore, at 'late, steady times', Ėloss approximately balanced Ėoxid and the system approached a steady-state condition.Towards steady-state conditions, the net stored energy approached a maximum, constant value (i.e., E net , seen as the red line in Fig. 2c).Note that nearly all stored energy accumulates in the cooling zone in applied smouldering systems [29,65].While Fig. 2 demonstrates this transition to steady-state conditions, no study has yet generalized the conditions that govern the characteristic energy evolution of smouldering systems.
This study presents a framework to describe the characteristic evolution of the global energy balance common to applied smouldering systems.This was achieved by deriving a new analytical model of the steady-state temperature profile throughout the cooling zonewhich accounts for nearly all stored energyto approximately define the maximum amount of stored energy.This framework was applied to a validated numerical model and experiments under various conditions (e. g., reactor width, fuel type, fuel concentration, insulation quality, and air flux).The results were made non-dimensional and shown to collapse onto a unified curve, which demonstrates the characteristic global energy balance evolution.Altogether, this analysis elucidates the key influences on the global energy balance, which can be manipulated by engineers in designing smouldering systems for environmentally beneficial purposes.

Expressing the energy balance with non-dimensional terms
To develop a characteristic energy balance curve, the transient net stored energy (E net ) and time (t) need to be expressed nondimensionally.The energy needs to be non-dimensional (E net ) to the maximum stored energy at steady-state conditions, i.e., E net,f when Ėoxid ≈ Ėloss in Eq. ( 1) [29,48,60], and the time needs to be non-dimensional (τ) to the approximate time required to reach this condition (t f ).These terms also need to account for the ignition effects and can be described as: where t ig is the ignition time (i.e., identified as when the heater was turned off), l ig is the heated length upon ignition, l cool,f is the steady-state cooling length, v oxid is the smouldering velocity, v cool is the cooling velocity [67], and E net,ig is the stored energy accumulated prior to ignition.
In convective and radiative ignition systems, l ig equals the smouldering distance upon ignition; in conductive ignition systems, l ig equals the smouldering distance plus the hot sand below the heater element upon ignition (which was approximately 0.1 m in LAB S1 and B1).In all cases, the l ig and l cool,f terms represent approximate lengths of hot regions.
Because nearly all the stored energy in applied smouldering systems using inert media accumulate in the reaction-less cooling zone [29,59], this zone was targeted to estimate the terms in Eqs. ( 3) and (4).

Estimating E net,f from a steady-state profile
To estimate E net,f , a new analytical solution was developed to describe this temperature field throughout the cooling zone based on the system and smouldering conditions.That is, the temperatures throughout the reaction and preheating zones were not modelled.This formulation followed the perturbation method from Kuznetsov [67,68] for local thermal non-equilibrium (LTNE) heat transfer between the gas (air) and solid (sand) phases and fixed the frame of reference to the smouldering front (i.e., common to smouldering analyses [58,61,62,69]).
Equations ( 5) and ( 6), were then transformed by fixing the frame of reference on the smouldering front, as the cooling zone was assumed to be fixed relative to the smouldering front, i.e., x = x − v oxid t [58].The transient term is then dropped to focus on steady-state cooling zone characteristics upstream of the smouldering front: Equations ( 8) and ( 9) were then made non-dimensional following simplifications and non-dimensional descriptions: where ṁ′′ bulk is a negative value under the conditions investigated, as C p ṁ'' s > C p ṁ'' g .Furthermore, δ was defined as very small parameter to permit the perturbation analysis [67,68]; therefore, the perturbation analysis is invalid at small reactor radii, i.e., ≲O (0.01 m).Equations ( 8) and ( 9) then become: where The following boundary conditions were then applied to Eqs. ( 8) and ( 9): ∂T g ∂r ∂T g ∂x These boundary conditions and key assumptions are illustrated in Fig. 3. Equation ( 14) embeds the assumption that smouldering proceeded with constant intensity across the reactor radius.While constant peak temperatures are often observed in robust smouldering systems [22,23,26,27,33,57], the peak temperatures are generally cooler near the wall due to radial heat losses [58,[74][75][76].Furthermore, near-wall extinction (i.e., due to non-uniform reactions [51,58,64]) and airflow divergences (i.e., due to non-uniform air flux [58,74,75]) also lead to complicated dynamics.As non-uniform air flux and reactions were not critical to describe the evolution of E net , they were not pursued.However, they are worthwhile to explore in multidimensional numerical models (e.g.Refs.[12,13,49,74,77]), as Rashwan et al. [26], demonstrated that they can lead to operational challenges, e.g., lower processing rates, peak temperatures, and heating rates.
Equation (15) describes the convective boundary at the edge of the insulation, where the temperature was assumed to decrease linearly across the thickness of the wrapped insulation.These assumptions are encapsulated in the modified heat transfer coefficient, , where r oins is the outer radius of the insulation [29].As the heat transfer coefficient (h ins ) at the outer radius of the insulation and the effective bulk conductivity of the insulation (k bulkins ) were unknown, H was varied as a bulk parameter in the sensitivity analysis in Section 3 (following [29]).Equation ( 16) assumes that the temperature did not vary axially far downstream of the smouldering front, which matches observations in long smouldering systems [48].The boundary conditions in Eqs.14-16 were then made nondimensional using Eq. ( 10): where R = r o /r o = 1 and α I = Hr o .Note that the sign of ∞ changed from Eqs. ( 16)-( 19) because of the sign change in the axial transformation from x to ξ x in Eq. (10).
With the boundary conditions in Eqs.17-19, the analytical solution to Eq. ( 12) was achieved following the methodologies outlined in Refs.[70,78] (where additional insight can be found in Refs.[79,80]).The solution is therefore: where J 1 (γ m ) = α I RJ 0 (γ m ).The difference in phase temperatures, Δθ(ξ x , ξ r ), was then solved using Eqs.( 13) and ( 20): Altogether, Eqs.20 and 21 are applicable to describe steady-state cooling zone temperatures across many applied smouldering systems with various fuels, fuel concentrations, air fluxes, reactor sizes, and insulation qualities.However, these systems need to facilitate sufficiently long and wide, reactionless cooling zones trailing a steady smouldering front (i.e., with constant peak temperatures and velocities) in systems with constant, uniform air fluxes and minimal LTNE.Therefore, Eqs.20 and 21 cannot describe relatively short cooling zones in small radii reactors, or in systems that are highly affected by nonuniform air fluxes, reactions within the cooling zone, large temperature differences between phases, or any other dynamics beyond convective LTNE heat transfer.

Estimating l cool,f from a global energy balance
While Eqs. 20 and 21 provide descriptions of the gas and solid phase temperatures throughout the cooling-zone at steady-state conditions, Eqs.22-24 provide a simplified description of the approximate final cooling zone length (l cool,f ) at the same conditions (i.e., when Ėoxid ≈ Ėloss ).This approximation was presented in Rashwan et al. [29], and is summarized below: T.L. Rashwan et al.where ΔH oxid eff is the effective heat of smouldering (which accounts for pyrolysis and other endothermic processeslike water phase change), m fuel /m s is the mass ratio of fuel to sand, x c is the characteristic heat transfer length from the convective-diffusive equation (i.e., x c = k bulk /ρ g C pg u g ) [67], and α II is the modified heat transfer coefficient used in Ref. [29] (i.e., α II = Hx c = α I x c /r o ).

Experimental and numerical results
The analyses from Sections 2.1-2.3 were applied to the data summarized in Table 1.DRUM R0-2, LAB R1-2, and LAB2 R0a-R0b were robust experiments that used granular activated carbon (GAC) in sand; LAB S1 used sewage sludge in sand; and LAB B1 simulated bitumen in sand using a 1D validated numerical model.The LAB experiments and LAB B1 used and simulated, respectively, reactors with 0.080 m radiusand the LAB2 and DRUM experiments used reactors with 0.054 and 0.300 m radii, respectively.See relevant schematics in the Supplementary Materials, Fig. S1.
Because all experiments and LAB B1 in Table 1 used and modelled, respectively, similar coarse grain sand, the same expressions for the air and sand thermophysical properties were used in the analyses in Sections 2.1-2.3.All effective thermophysical parameters were estimated over T amb to T peak , e.g., ρC p eff = ∫ T peak T amb ρC p (T)dT/(T peak − T amb ) [26,29,58].
See Table 2 for all parameters and references [27,56,60] for thorough descriptions of the parameter errors.
As noted in Table 1, all results have been detailed elsewhere, except LAB S1, which was an experiment with 72% MC sewage sludge mixed with coarse grain sand in a 1.135 m tall fuel bed.This experiment followed established experimental methods for conductive ignition [23,26], and was instrumented with centreline thermocouples (spaced 0.02 m apart from 0.01 m below the heater to the top of the fuel bed) and continuous emissions analysis above the top of the pack (for O 2 , CO 2 , and CO).The same equipment detailed in Refs.[23,66] was used in this experiment (e.g., thermocouples, datalogger, insulation, conductive heater, coarse grain sand, source of sewage sludge).This experiment facilitated self-sustaining smouldering throughout most of the reactor but weakened and developed a large crust formation (further detailed below in Section 3).LAB S1 eventually quenched after smouldering 0.85 m.It is hypothesized that this quenching was due to water condensation and accumulation ahead of the smouldering front; this phenomena was observed in smouldering experiments using surrogate faeces under similar conditions [22].Therefore, this experiment was analyzed prior to quenching when it achieved an approximately steady-state centreline temperature profile in the cooling zone.Note that these steady-state conditions are seldom observed in applied smouldering systems because they are often too short to allow cooling zone to fully developas it can grow multiple metres in robust systems based on common ranges of the terms in Eq. ( 24) [29].All other experiments aside from LAB S1 in Table 1 reflect these typical smouldering conditions and were unsuitable for steady-state analyses.Therefore, LAB S1 provides unique and valuable experimental data to understand steady-state conditions within the cooling zone.

Details on ΔH oxid eff estimations
The ΔH oxid eff values from DRUM R0-R2, LAB R1-R2, and LAB R0a-b (which used GAC) were estimated from the chemical equation: [58,60,74,75].That is, the combustible fraction of GAC was assumed to be entirely carbon and degrade into CO and CO 2 , and corrected for the GAC ash and moisture contents (see additional details and evidence in support of this procedure [27,60]).Therefore, 25 MJ kg − 1 was used as ΔH oxid eff for DRUM R0-R2, LAB R1-R2, and LAB R0a-R0b based on the measured CO/CO 2 fractions [27,58,73].20 MJ kg − 1 was used as the ΔH oxid eff for LAB B1, which accounted for the simulated exothermic and endothermic reactions from smouldering bitumen [48].As the sewage sludge used in LAB S1 is a complex fuel that is influenced by many unresolved mass transfer processes and reactions, ΔH oxid eff could not be simply estimated from CO/CO 2 signals or from known heats of oxidation and pyrolysis.Instead, the heat of reaction on a dry mass basis (ΔH oxid,dry ) for sewage sludge was estimated from Eq. ( 25): (1 − MC crit )ΔH oxid,dry = MC crit ΔH vap (25) where MC crit is the critical moisture content that leads to quenching (i.e., observed experimentally as approximately 80% from Rashwan et al. [23]) and ΔH vap is the heat of water vaporization from 25 • C (2.45 MJ kg − 1 ).Therefore, upon correcting for the MC in LAB S1 (i.e., 72%), the ΔH oxid eff for LAB S1 was estimated as 0.98 MJ kg − 1 following: (26) Note that the energy balances in Eqs. 25 and 26 are highly simplified and a more rigorous analysisif availablewould be beneficial to estimate ΔH oxid eff .However, the authors are not aware of a more reliable analysis with the data available from LAB S1.ΔH oxid eff is notoriously difficult to estimate in smouldering systems [27,46].Therefore, in absence of an alternative analysis, Eqs. 25 and 26 provide a useful estimate of the energy released from LAB S1 to understand the steady-state conditions in that experiment.

Estimating E net /A profiles from data
While the analyses outlined in Sections 2.1-2.3 provide a framework to fully describe the energy balance evolutions in smouldering systems, many smouldering studies only instrument or simulate the centreline (r = 0).Therefore, to compare with all data available in Table 1, only the centreline temperature data was used to define the stored energy per cross-sectional area along the centreline, E net /A [27,29]: where x 0 and x f bound the fuel bed.The theoretical maximum stored a Average temperature dependent values are examples over the centreline T amb to T peak measured in DRUM R2.
T.L. Rashwan et al. energy in the cooling zone, E net,f /A, was estimated as: where θ s = θ g + δΔθ is the solid phase temperature modelled along the system centreline from Eqs. 20 and 21 and 2l cool,f from Eq. ( 24) was used as a practical boundary for Eq. ( 28) to delineate the end of the cooling zone; most of the cooling occurs between 0 ≥ x ≥ − l cool,f but the integration in Eq. ( 28) was extended − 2l cool,f to account for trailing temperatures (see an examples of l cool,f alongside analytical profiles in Fig. 5).
Following the results from Rashwan et al. [29], α II = 0.5 was used to estimate E net,f /A and t f values for DRUM R0-R2 and LAB R1-R2.Moreover, because the LAB2 experiments from Zanoni et al. [73], used better quality insulation than DRUM R0-R2 and LAB R1-R2, α II = 0.15 was also used to estimate E net,f /A and t f values for LAB2 R0a-b [65].The respective reactor radii were used as the r o terms for all analyses, except for LAB S1.LAB S1 exhibited a restricted effective radius due to a large crust formation; therefore, α II = 100 and r o = 0.055 m were used to estimate E net,f /A and t f from that experiment (see justification below in Section 3).The E net /A measured when the heater was turned off (t = t ig ) was used as E net,ig /A in all analyses.To make the results from LAB B1 non-dimensional: the final E net /A measurement was used as E net,f / A (i. e., from the maximum E net value in Fig. 2c) and t f was estimated with Eq.
Though E net,f /A accounted for most of the stored energy in the systems (i.e., in the cooling zone), it did not account for the small amount of energy stored in the reaction and preheating zones (Fig. 1), which were captured in E net /A valuessee Eq. (27).Therefore, at near steady-state conditions, E net /A is larger than E net,f /A.Altogether, many simplifying assumptions were used to approximate E net,f /A and t f , which are therefore imperfect but useful.Regardless of the limitations, the key energy balance trends can be identified by making E net / A profiles nondimensional following: The non-dimensional energy profile (E net /A) can then be evaluated from various data sources and compared together.

Results and discussion
Fig. 4 compares centreline temperature spatial profiles from three example experiments (DRUM R2, LAB R2, and LAB S1 in Fig. 4b, c, and  d, respectively) alongside a theoretical adiabatic profile in Fig. 4a.The temperatures were made non-dimensional following Eq.( 10) and the distances were normalized to the cooling zone length (l cool ) at each snapshot in time.The data from DRUM R2 and LAB R2 were gathered near the end of smouldering, i.e., when the smouldering front was near the end of the fuel bed (x f ) respective to each system, and from LAB S1 after steady-state conditions were reached (τ > 1).Note that the temperature markers in Fig. 4d ahead of smouldering (i.e., righthand side) are not at ambient values because of evaporation and recondensation of the wet sewage sludge [54].The sizes of the smouldering and cooling velocity vectors (v oxid and v cool , respectively) indicate their relative speeds.
The E net /A values in: (i) Fig. 4a is a theoretical minimum in an adiabatic system; (ii) Fig. 4b-c were the maximum values measured; and (iii) Fig. 4d is a theoretical maximum in a steady-state system, which was close to the E net /A values measured (see LAB S1's E net /A profile in Fig. 7).Fig. 4 also illustrates the centreline temperature profile corresponding to changes in E net /A.Fig. 4a shows a theoretical condition where all energy generated is stored in the system.In this case, E net /A = 0 because E net,f /A tends to infinity in an adiabatic system, as the cooling length grows unboundedlysee Eq. ( 29).Fig. 4b shows the temperatures from DRUM R2, which was an experiment that smouldered GAC in a large reactor (0.3 m radius) with minimal heat losses relative to the energy released after propagating 0.8 m with E net /A ≈ 0.1.Fig. 4c shows the temperatures from LAB R2, which was an experiment that followed the same conditions as DRUM R2 but in a small reactor (0.08 m radius), after propagating 0.5 m, where E net /A ≈ 0.4.In other words, the increase in E net /A from Fig. 4b to c indicate the increased influence of heat losses due to the decreasing reactor radii [27,29].Fig. 4d shows a key limiting case from LAB S1 after propagating 0.5 m with E net /A ≈ 1.This high E net /A value indicates that the stored energy was near the system's maximum value when Ėloss and Ėoxid roughly balancesee Eq. (1).Note that the l cool (t) values on Fig. 4a-c were transient and grew throughout propagation, as v oxid > v cool [29,60].In these frames, v cool was controlled by the axial heat transfer, i.e., v cool = ρ g C pg u g /(ρC p ) bulk [67].
On the other hand, the cooling zone length in the steady-state condition Fig. 4. (a-d) Conceptual models (solid lines) of the centreline spatial temperature profiles transitioning from (a) adiabatic (E net /A = 0) to (d) approximately steadystate conditions (E net /A = 1) when the energy released rate from smouldering approximately equals the heat loss rate from the system, i.e., Ėoxid ≈ Ėloss .The overlain data markers are from snapshots from select experiments that demonstrate this transition.The black, grey, and white markers in (b), (c), and (d), respectively, are the centreline temperature profiles from DRUM R2, LAB R2, and LAB S1, respectively (see Table 1).The net stored energy along the centreline (E net /A), smouldering velocities (v oxid ), cooling velocities (v cool ), and cooling zone lengths (l cool ) are labelled for further discussion in the text; note that the maximum final cooling zone length and cooling velocity are l cool,f and v cool,f , respectively.
T.L. Rashwan et al. in Fig. 4d was constant at a maximum value, i.e., l cool,f , which can be estimated from Eq. ( 24).Moreover, the cooling velocity in this condition (v cool,f ) was also a maximum value and controlled by perimeter heat losses, where v cool,f ≈ v oxid .Altogether, Fig. 4 illustrates the characteristic evolution of spatial temperature profiles in applied smouldering systems, from adiabatic (E net /A = 0) to steady-state (E net /A ≈ 1) conditions.
Fig. 5a-b compares temperatures measured from LAB S1 with those modelled with Eqs.20-24 to build confidence in the model formulation, where solid and dashed lines indicated air and sand temperatures, respectively.The temperature measurements (circles) were taken when the smouldering front travelled 0.71-0.81m, as identified via the peak temperatures over this period, which averaged 496 • C.This timeframe was chosen so to show a wide spread of the trailing cooling zone temperatures that were well-within LAB S1's late-time period (i.e., τ > 1), which occurred after the smouldering front travelled approximately 0.26 m.The smouldering front was anchored at x = 0 m in Fig. 5a-b to compare with Eqs.20 and 21.Note that the steady-state cooling lengths (l cool,f ) from Eq. ( 24) were compared alongside the modelled temperatures at corresponding model conditions (i.e., horizontal lines below the temperatures).Sensitivity analyses were completed on the insulation quality (α II ) and the outer radius (r 0 ) to explore the influence of these key system parameters, and better understand the consequences of the crust formation in LAB S1.Fig. 6 is a top-down photo of the cross-section from LAB S1 upon excavation.
Fig. 5a shows the small influence of α II from 1 to 100 on the temperatures with r 0 = 0.08 m.This insensitivity indicates that when the insulation is poor, it effectively leads to a constant temperature condition at the outer radius, i.e., heat transfer out of the insulation is fast and the outer-radius remains near T amb (see Supplementary Materials, Fig. S2).This range of poor insulation (i.e., high α II values) was explored to approximate the effects of a crust formed near the wall from LAB S1.As shown in Fig. 6, LAB S1 exhibited a large, unburned ring of crust.This crust drew large heat losses that were not explicitly modelled, only approximated with high α II values.
While the model results in Fig. 5a do not match with experimental measurements, Fig. 5b shows much better matching by decreasing r 0 with α II = 100.Though LAB S1 was performed in a 0.08 m radius reactor, Fig. 6 shows that the unburned edge grew non-uniformly inward, which effectively reduced the radius for smouldering and energy storage in the cooling zone.From the sensitivity analysis in Fig. 5b, the modelled temperatures with an effective radius near 0.05-0.06m best captured the observed temperatures, which approximately match the observed crust growth in Fig. 6.Therefore, r o = 0.055 m was used to estimate t f and E net,f /A from LAB S1.In Fig. 5a-b, the overlain sand and   air temperatures show comparable differences as previous perturbation analyses with similar analytical models [26,29,58].In Fig. 5a-b, these temperature differences increase as α II increases and r o decreases, which indicates that phase temperature differences dampen in long cooling zones with gradual temperature gradients.Moreover, the elevated solid temperatures near x = 0 m at small radii in Fig. 5b were due to solution instabilities, e.g., most noticeable with r o = 0.04 m.The perturbation model is least reliable at small radii (as discussed in Section 2.2).Moreover, the l cool,f values in Fig. 5a-b approximately captured the cooling zone lengths of the modelled temperature profiles.Altogether, Fig. 5 illustrates that the model assumptions embedded Eqs.20-24 provide valuable insight into applied smouldering systems' cooling zones at steady-state conditions.Fig. 7 illustrates the E net /A profiles during smouldering from all models and experiments in Table 1 over non-dimensional time (τ) via the methods summarized in Section 2. Note that no experimental curve perfectly approaches E net /A→1 because many simplifying assumptions were needed to harmonize the data (see details in Section 2).Moreover, the experimental curves exhibit false oscillations, which were primarily due to coarse thermocouple placements [27].As an exception, the LAB B1 exhibits a smooth curve and does reach E net /A = 1 at τ = 2.4, as the maximum E net /A value modelled was used as E net,f /A (see details in Section 2.6).Most importantly, all E net /A profiles generally collapse onto each other; therefore, Fig. 7 illustrates that the assumptions used throughout the analyses in Section 2 do successfully clarify the non-dimensional, characteristic energy balance evolution of applied smouldering systems.
As discussed throughout this manuscript, the experiments and models exhibited many differences; however, the implications of these differences on their global energy balances can be distilled in the extent of their E net /A profiles.That is, DRUM R0-R2 show relatively short profiles in Fig. 7 that progress until E net /A ≈ 0.1, τ ≈ 0.03− 0.04 (see the inset in Fig. 7b for clarity).These short profiles illustrate that DRUM R0-R2 operated only in early durations of the characteristic energy curve.In other words, the global energy balances from Eq. (1) in DRUM R0-R2 were far from steady-state conditions and these experiments accumulated large amounts of released reaction energy throughout their operation.Conversely, the other experiments and simulations in Fig. 7 exhibited longer profiles than the DRUM experiments.For example, the profiles from LAB R1-R2 and LAB2 R0a-R0b were longer as they used smaller radii (i.e., DRUM: 0.300 m, LAB: 0.080 m, and LAB2: 0.054 m) at similar smouldering conditions (e.g., fuel type, concentration, and applied air fluxsee Table 1).The decreasing reactor radii promoted greater relative heat losses, which drove the global energy balances in LAB R1-R2 and LAB2 R0a-R0b closer toward steady-state conditions [27,29].The simulation LAB B1 and experiment LAB S1 both exhibited the longest profiles on Fig. 7, which approached E net /A→ 1. LAB B1 and S1 demonstrate global energy balance evolutions towards steady-state conditions when the rate of heat losses ( Ėloss ) approximately balance the rate of energy released from exothermic smouldering reactions ( Ėoxid ).While Zanoni et al. [48], demonstrated this evolution numerically under specific conditions in Fig. 2, this study identifies the key parameters that govern this transition generally.Altogether, Fig. 7 demonstrates the characteristic global energy balance evolution in applied smouldering systems, which can be manipulated by engineers to better design future systems.

Conclusions
Under a broad range of experimental and modelling conditions, the global energy balance evolutions in applied smouldering systems have been examined and a unified, characteristic profile has been elucidated.This characteristic profile was uncovered through using a description of the global energy balance at steady-state conditions, i.e., when the rate of energy released from smouldering approximately balanced perimeter heat losses.This analysis exploited the fact that most stored energy in applied smouldering systems accumulate within the cooling zone.Therefore, a new analytical description of cooling zone steady-state temperature distribution was used to estimate the maximum stored energy, and the cooling and smouldering velocities were used to approximate the time needed to reach these conditions.
Altogether, this study provides novel insight into the key drivers that govern the characteristic energy distribution and evolution in applied smouldering systems, which can be manipulated by engineers and researchers to design improved applied smouldering systems.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.(a) Conceptual model of the temperature distribution through the centreline of a smouldering system (solid line), overlain with temperature data (markers) from the end of propagation through a 0.9 m tall mixture of granular activated carbon and coarse grain sand (i.e., the end profile of DRUM R2, detailed in Table1).The key zones: Inert Heating Zone, Reaction Zone, and Cooling Zone are highlighted for further discussion and are approximately bounded between the velocities of the smouldering front (v oxid ) and the cooling front (v cool ), where v oxid > v cool .(b) A picture of the oil-drum sized reactor (DRUM) used to conduct the experiment plotted in (a).

Fig. 2 .
Fig. 2. LAB B1 numerical modelling results from Ref. [48] that show: (a) temperatures, (b) the terms global energy balance terms from the righthand side of Eq. (1), and (c) the balance of Eq. (1) terms (dE net /dt) and the net stored energy (E net ).The profiles in (a-c) are normalized over dimensionless time (DT)[66], where DT = 0 notes ignition and DT = 1 notes the end of propagation (i.e., when the smouldering front reached the end of the fuel bed, which is marked in (a-c) with a dashed blue vertical line).The horizontal dashed grey line in (c) notes Ėnet = 0 and the dashed vertical black line in (a-c) approximately shows when Ėloss (t) = 0.9 Ėoxid (t) (i.e., when the front propagated around 1.6 m just after DT = 0.5).The black vertical line delineates the 'early, transient times' from the 'late, steady times'.This figure was adapted from Refs.[29,48], with permission of Elsevier.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 3 .
Fig. 3. (a) Schematic of the smouldering system where (i), (ii), and (iii) show the Inert Heating Zone, Reaction Zone and Cooling Zone, respectively; and (b) a simplified conceptual model with the boundary conditions used to solve the temperature profile within the cooling zone.This figure was adapted from a similar schematic in Rashwan et al.,[29], with permission of Elsevier.

Fig. 5 .
Fig. 5. (a) and (b) cooling zone temperature data from LAB S1 (circles) compared to the modelled sand (dashed lines) and air (solid lines) temperatures from Eqs. 20 and 21 under various model assumptions.origin in (a) and (b) is anchored to the peak temperature observed as the smouldering front, and the negative coordinates indicate the cooling zone behind the smouldering front.The straight lines underneath the analytical curves indicated the maximum cooling zone lengths, i. e., l cool,f from Eq. (24), at the same conditions as the analytical curves.The sensitivities to insulation quality (α II ) and reactor radius (r o ) are overlain to (a) and (b), respectively.Note that r o = 0.08 m and α II = 100 in (a) and (b), respectively, and the α II = 10 and α II = 100 curves are nearly overlain in (a).

Fig. 6 .
Fig. 6.A photo of the crust from LAB S1 with dimensions labelled.

Fig. 7 .
Fig. 7. (a) The non-dimensional energy per area and time results from all data in Table 1.Only every tenth data point from the experimental results (solid lines) were plotted for clarity.All data points from the simulation LAB B1 (dashed line) were plotted.(b) The inset isolates the profiles from DRUM R0-R2.

Table 1
[23,57,66,81,82]ditions and key smouldering front results.Experimental Conditions' errors represent conservative estimates of experimental equipment error.Smouldering Front Results' errors are conservative estimates of experimental variability and are normalized standard deviations from three DRUM and LAB repeat experiments from smouldering wastewater sewage sludge.These errors align well with similar experimental studies[23,57,66,81,82].The LAB errors are representative of the LAB2 systems, which used the same equipment except a different sized reactor.
[27] smouldering front metrics capture smouldering conditions sufficiently far from ignition and end effects[27].b At standard temperature and pressure (21.1 • C at 1 atm).c The percentages represent volume fractions and fr CO = CO% /(CO% + CO 2 %).d CO measurement throughout propagation exceeded DRUM emissions analyzer calibration range, i.e., 0-0.3 vol%.e CO production was not simulated.

Table 2
Temperature-dependent model input parameters.