Analytical and numerical insights into wildfire dynamics: Exploring the advection–diffusion–reaction model

Understanding the dynamics of wildfire is crucial for developing management and intervention strategies. Mathematical and computational models can be used to improve our understanding of wildfire processes and dynamics. This paper presents a systematic study of a widely used advection–diffusion– reaction wildfire model with non-linear coupling. The importance of single mechanisms is discovered by analysing hierarchical sub-models. Numerical simulations provide further insight into the dynamics. As a result, the influence of wind and model parameters such as the bulk density or the heating value on the wildfire propagation speed and the remaining biomass after the burn are assessed. Linearisation techniques for a reduced model provide surprisingly good estimates for the propagation speed in the full model


Introduction
Wildfires are a crucial part of the Earth system that have shaped ecosystems since the appearance of terrestrial plants [6].While many ecosystems depend on a certain wildfire regime to sustain themselves [17], wildfires pose a significant threat to anthropogenic biomes, destroying property and, in the worst case, leading to casualties.In addition, large and uncontrolled wildfires are often accompanied by air pollution and significant change of the hydrological function of the impacted area that may harm water security [22].At the time of writing, a large wildfire in Nova Scotia, Canada, is making headlines as it burned down a record area of about ten million hectares.A huge wildfire is ravaging the Greek island of Corfu, another wildfire at the coast of Algier claimed 34 lives, and yet another wildfire is approaching the city of Palermo, Italy.In previous years, Californian [21,22] as well as Australian wildfires [10] wreaked havoc on the regions.The Mediterranean region that has been historically prone to wildfires, is scorched by early wildfires, see, for example [32,34].South East Asia is reportedly bracing for increasing wildfires, see for example [26].The African Sahel region is more frequently subject to wildfires and agricultural fires that may in return trigger prolonged droughts [19].For a current global review of wildfire dynamics, see the excellent discussion in [17].
In general, the anthropogenic climate change continues to create favorable conditions for wildfires and both frequency and intensity of wildfires are expected to increase [12,28,33].Previously unaffected regions such as Northern Europe have also become vulnerable to wildfires due to climate change [9].As wildfire risk increases, the demand to better understand and predict wildfire dynamics through mathematical models and their numerical solution to support management decisions is increasing.
Wildfire dynamics consist of two main processes: (i) ignition and (ii) propagation [17].Neither of these two processes are completely understood [17].In this context, mathematical and computational models can help to understand key dynamics of the fire-vegetation relationship [11].A comprehensive review of the modelling of wildfire propagation dynamics can be found in [29][30][31].In this work, we focus on a specific class of physically-based wildfire propagation models, the advection-diffusion-reaction equation-based (ADR) model.This ADR model has been initially proposed in [2] for wildfire propagation in an idealised two-dimensional (2D) forest.It was numerically explored in [7,8].Mathematical properties-like the existence of travelling wave solutions sensu [16]-were studied for simplified models in [4].
Variations and extensions of this ADR wildfire model exist in the literature.For example, in [5], a threedimensional (3D) form of the ADR wildfire model is presented that allows to distinguish between crown and surface forest fires.In [3], the model was extended with convection above the fire by means of coupling the ADR wildfire model with an atmospheric flow model.Another variation of the ADR wildfire model is discussed in [14], which explicitly accounts for the underlying chemical reactions of combustion and spread in the framework of an ADR.In [23], an ensemble Kalman filter technique is used to account for measured temperatures in running ADR-based wildfire simulations.
While the ADR wildfire model is widely used and adapted, the individual model components and their influence on the dynamics have not yet been systematically studied.In this paper, we formulate sub-models of the ADR wildfire model, each including combinations of the mechanisms advection, diffusion, reaction and ambient cooling.We study the influence of these hierarchical sub-models on the full ADR wildfire model's emerging non-linear behaviour using both analytical and numerical methods.The rest of this paper is structured as follows: In Sec. 2, the ADR wildfire model is derived using physics-based arguments and the numerical solver used to solve the resulting system of equations is presented.The sub-models are analysed in Sec. 3. The interplay of all mechanisms is studied in Sec. 4, first for a one-dimensional domain and then for a two-dimensional more realistic setting.Finally, conclusions are drawn in Sec. 5.
2 The wildfire propagation model

Advection-diffusion-reaction equation-based wildfire model
The model presentation in the literature, for example [2,7,8], is focused on mathematical aspects.Here, we provide a physics-based description of the model.The model is characterised by the biomass fraction Y (x, t) and the temperature T (x, t), defined inside a two-dimensional spatial domain Ω ⊂ R 2 , with x = (x, y) ∈ Ω being the Cartesian coordinates in space and t > 0 being time, see Fig. 1.The temperature is a measure for the internal energy and we will derive the model by energy considerations.Bulk parameters are used to characterise the properties of the biomass and air mixture, which are considered as a continuum.
Our starting point is a fixed control volume defined as V = S × [0, ℓ] inside a forest ecosystem.Here, S ⊂ Ω is the horizontal extent of the forest and l is the forest layer thickness or depth in the vertical direction, associated with the tree height in the forest.The boundary of the control volume is denoted by ∂V .The equation for the conservation of internal energy inside V is where e is the specific internal energy per unit mass, ρ is the density defined below, and Qcond , Qconv , Qrad and Qreac are the conduction, convection, radiation and reaction heat fluxes, respectively.A conceptual sketch of these fluxes is shown in Fig. 1.Note that compressibility effects and viscous dissipation in the gaseous phase are neglected in the energy balance.The advection velocity, v, is related to the wind velocity.
We assume that all modelled quantities are homogenised over the forest layer thickness l and thus, represent depth-averaged quantities on a two-dimensional horizontal plane defined by S. Consequently, heat fluxes and wind velocity only have horizontal components.The fuel mass fraction Y relates to the density ρ as where m f is the fuel mass and m f0 is the maximum fuel mass at t = 0. Here, we assume that the total mass of air in V is negligible such that the total mass m t = m a + m f ≈ m f , and therefore ρ ≈ ρ f .Assuming a constant density ρ 0 and specific heat c in V , the internal energy can be expressed as e = cT .Therefore, the temperature is a quantity for the internal energy.This allows us to rewrite Eq. ( 1) as The conduction heat flux Qcond is given by Fourier's law as where k is the thermal conductivity.The convection heat flux is given by Newton's cooling law with h 0 being the convection coefficient and T ∞ being the temperature of the air over the surface.The radiation heat flux can be written in the form of a non-linear conduction heat flux as provided that the optical path length for radiation δ is smaller than the characteristic length of the control volume.Finally, the reaction heat flux is given by where H is the combustion heat per unit mass of fuel and is the fuel mass disappearance rate per unit volume, with being the rate of variation of the mass fraction of fuel, given by the Arrhenius law where T ac = E A R with the activation energy E A , the universal gas constant R, a pre-exponential factor A and the activation function s(T ) given by that triggers the reaction when the temperature is higher than T pc .By combining the previous definitions and inserting them into Eq.( 7), we obtain Using the Green-Gauß divergence theorem and considering that the flow is incompressible, so ∇ • v = 0, we can rewrite Eq. ( 3) as with h = h 0 /l, which allows to obtain Finally, the complete wildfire propagation model is written as where we define k t (T ) = 4σδT 3 + k and Ψ(T ) = s(T )A exp − Tac T .The temperature T = T (x, t) > 0 and the biomass Y = Y (x, t) ∈ [0, 1] are the model variables.The model parameters are summarised in Tab. 1 with values taken from the literature when available.Parameters that could not be assigned a value range based on the literature are estimated, often through heuristic arguments.In general, many of these parameters may vary in space due to spatial heterogeneity in the environment.Examples for potentially spatially heterogeneous parameters are the bulk density ρ 0 and the specific heat c.In this work, we will not consider spatial variations of this kind.To complete the mathematical model, Eq. ( 14) must be subjected to appropriate initial and boundary conditions, IC and BC, respectively.Many choices exist, for example IC: The zero-flux boundary condition defined in Eq. ( 15) could be applied if the domain is chosen such that the wildfire does not propagate outside of it.

Numerical solver
The ADR wildfire model outlined in Eq. ( 14) and Eq. ( 15) is discretised using a high order finite volume scheme.The weighted essentially non-oscillatory (WENO) method [20] with up to 7-th order of accuracy is used to compute the convective fluxes.The WENO method provides stable reconstructions of cell variables without Gibbs oscillations in the presence of sharp gradients or discontinuities in the solution, which may be the case for the problems herein considered.Second order central differences are used for the diffusive fluxes.This combination proved to be suitable for other applications in [24,25].Different Runge-Kutta integrators of varying order are used in combination with Strang splitting to step forward in time, similar to the numerical solver in [7].The time integration is as follows: First, the reaction Table 1: Summary of model parameters and their range used in this study with references when applicable.
terms in Eq. ( 14) are evaluated to update the model variables for half a time step using an implicit second order Runge-Kutta (RK2) method.Then, the advection and diffusion terms in the partial differential equation are evaluated to evolve the system for one full time step using a third order Strong Stability Preserving Runge-Kutta (SSP-RK3) method.Finally, the reaction terms are evaluated again for half a time step using the RK2 method.The validation of the numerical solvers for a spatially integrated model and for the spatially explicit ADR wildfire model are presented in Appendix A.1 and A.2, respectively.The need of very high order of accuracy in space for the discretisation of the convective fluxes is also motivated in Appendix A.2.

Analysis of hierarchical sub-models
The ADR wildfire model in Eq. ( 14) includes the interacting mechanisms advection, diffusion, and reaction that lead to highly non-linear behaviour.See, for example [27] for a discussion of the manifold behaviour of reaction-diffusion equations.In order to better understand the effects of these individual mechanisms, we formulate sub-models of system (14) with fewer mechanisms.All models together form a model family and every model will provide insights into some characteristic behaviour of the system.In the following, for every sub-model, we summarise relevant results from the literature, prove analytical results, and carry out numerical simulations.We interpret the analytical and numerical results in the light of wildfire management applications and specify open questions.

Dynamical system analysis for the spatially lumped sub-model
We start by analysing the model in a spatially integrated or "lumped" version, neglecting any spatial mechanisms.The ADR wildfire model in Eq. ( 14) reduces to an ordinary differential equation model of the form with initial conditions T (0) = T 0 > T pc and Y (0) = Y 0 > 0. The combustion function Ψ(T ) is given by Ψ(T ) = s(T )A exp − Tac T with the discontinuous activation function s(T ) in Eq. ( 10) depending on the ignition temperature T pc .The discontinuity of s(T ) results in a discontinuity of the differential equation.
Proof.For T ≥ T pc the activation function s(T ) has the value 1.The right-hand side of Eq. ( 16) is continuous w.r.t.T and Y .Additionally, the right-hand side is differentiable w.r.t.T and Y , and the derivatives are bounded for the given domain.Consequently, the Picard-Lindelöf theorem guarantees the existence of a unique solution.
For T < T pc , the activation function is zero.For Y > 0 and T ∈ (T − ϵ, T + ϵ), ϵ > 0, the function f T (T, Y ) in T ′ = f T (T, Y ) is discontinuous.The existence of a unique solution is therefore not guaranteed and worth a proof.Lemma 2. The discontinuity of the functions f T and f Y for T = T pc does not affect the validity of the differential equations describing the combustion process.By interpreting the state (T pc , Y ) as new initial values for the differential equation the Picard-Lindelöf theorem guarantees the existence of a unique solution for the case T < T pc .
Consequently, the system given in Eq. ( 16) has a unique solution for given initial values When modelling the dynamics of a wildfire, the boundedness of the temperature by a lower and an upper bound is desired, where the latter may depend on the initial values.The system state Y describes a fuel mass density-compare Eq. ( 2)-and should be bounded therefore by Y ∈ [0, 1].We then prove that the solutions of system ( 16) are bounded.Assume The temperature T is bounded from above because Y is bounded by 1 and Eq. ( 18) gives This ensures a negative derivative for On the other hand, if T ∞ < T 0 < T pc , the temperature decreases with a lower bound T ∞ .These two cases neglect the activation of the combustion process.If, in the last case, T 0 > T pc > T ∞ yields, the combustion process starts and T ∞ again gives a lower bound for T .
Consequently, the solutions of System ( 16) are bounded, which is crucial for the interpretation of the solutions as a model for real world processes.Further, the switching behaviour of the combustion function leads to infinitely many stationary points, which differ in the biomass at ambient temperature.Proof.The change of biomass is only zero for either Y = 0 or s(T ) = 0.In the first case, for Y = 0, the temperature decreases until T = T ∞ and we find the stationary point (T ∞ , 0) T .The switching function is zero for T < T pc .Then, the derivative of the temperature is zero if T = T ∞ , which gives us infinitely many stationary states For gaining further understanding of the system dynamics, we visualize in Fig. 2 the phase space portrait of System ( 16) using the fixed parameters from Tab. 1.The heating and cooling regions in the phase space are indicated with red and blue colours.These regions are separated for T > T pc by the tipping line with T ′ = 0 given by The dashed gray tipping line in Fig. 2 corresponds to the maximal temperature for any trajectory.The continuous line of stationary points of the system is plotted with a green line in Fig. 2.
Remark 5.In the case where combustion takes place, so for initial values T 0 > T pc , only stationary points with Y ⋆ < Y tip (T ∞ ) are stable and occur as results of the systems dynamics.
We conclude the analysis of the spatially lumped model: Even if the equations contain noncontinuous terms, the existence of unique and bounded solutions is guaranteed.The set of stationary points is a continuous line with varying values for the biomass.Based on these analytical results, we now investigate the system numerically.
System ( 16) is solved for the initial conditions T 0 = 470 K and Y 0 = {0.2,0.4, 0.6, 0.8, 1.0} using the numerical solver described in Sec.2.2.The evolution in time of the temperature and biomass is depicted in Fig. 3.The solutions show the switching behaviour of the combustion function for T = T pc = 400 K resulting in an unsteady slope of the temperature decrease.The trajectories of these numerical solutions are also depicted in the phase space portrait in Fig. 2 using black continuous lines.Starting with the same initial temperature and various initial biomass values, the trajectories lead to comparable final values for the biomass, and the ambient temperature.The influence of the initial biomass on the remaining biomass after the combustion process is therefore small.This is more clearly visualised in Fig. 4, where the terminal biomass Y ⋆ is plotted for various combinations of initial temperature T 0 and initial biomass values Y 0 .Fig. 4 shows that a majority of the initial value combinations result in similar terminal biomass at the end of the simulation.
Fig. 5 shows the tipping line in Eq. ( 20) for different values of the relevant parameters.We consider as the base parameterisation the values in Tab. 1.In the plot, m = ρ 0 HA represents the denominator of h ρ0HA , where ρ 0 is the bulk density.A larger parameter m leads to a smaller upper bound for the remaining biomass.This can be interpreted as a larger bulk density leading to a higher temperature in the burning   process and therefore to a lower biomass in the end.A larger activation temperature T ac leads to a larger upper bound for the remaining biomass.This can be interpreted as a reduction of the availability starting a combustion process, perhaps linked to forest management decisions.
These effects on the tipping line are as well visible if we calculate the sensitivity of the tipping line on the parameters directly.The sensitivity of Y tip on ρ 0 is given by Enlarging the parameter ρ 0 therefore has a negative effect on the tipping line, shifting it downwards.The same is observed for the heating value H.In comparison, the influence of the activation temperature is positive, due to These expressions of the sensitivity therefore support the findings in Fig. 5.The tipping line in the phase space of System (16) gives an estimation for the maximum of the remaining biomass Y ⋆ after the combustion process.The minimum of Y ⋆ is an interesting quantity as well.As the system is non-linear and the dynamics are cut off by the switching function for T = T pc , a prediction of the steady states depending on the parameters and the initial values is challenging, even if System ( 16) is linearised.For example, a linearisation of the non-linear combustion term is where δ is chosen such that the linearised function is larger than 1 for T > T max .With this linearisation, the system reads and can be reformulated to The solution of the homogeneous equation is where the integration constant is a very large value, shifting the whole approximation above the domain of Y .This approach therefore is not giving useful bounds for the inhomogeneous and non-linear system.The analysis of the lumped System (16) shows the influence of some parameters on the tipping line and gives insight into the boundedness of solutions and the sensitivity of the model parameters.However, the minimum of Y ⋆ is further affected by diffusion and advection.In what follows, we analyse the effect of these mechanisms on the solution through sub-models of the ADR wildfire model.We start with combustion-free sub-models.

Combustion-free sub-models
Combustion-free sub-models of the ADR wildfire model neglect the combustion function Ψ in Eq. ( 14) to remove the discontinuity in the model.This allows to investigate paths to stationary states that result primarily from energy dissipation and the development of advection-driven travelling wave solutions.Here, we study two combustion-free sub-models: (i) a reaction-diffusion sub-model without combustion in Sec.3.2.1 and (ii) a combustion-free ADR sub-model in Sec.3.2.2.

Energy dissipation in a combustion-free reaction-diffusion sub-model
First, we consider a combustion-free reaction-diffusion model for x ∈ Ω ⊂ R d and constant diffusion.The reaction only describes the exchange with the ambient temperature and the word reaction is meant only in the context of reaction-diffusion equations.Without any combustion, the biomass is unchanged and Y ′ = 0. Further, we neglect advection, so v = 0.The system in Eq. ( 14) becomes with α = k t /(ρ 0 c) the thermal diffusivity and β = h/(ρ 0 c).We adapt the zero-flux boundary conditions from Eq. ( 15) and use initial conditions with T ≥ T ∞ point-wise and a constant Y (0, x) = Y 0 ∈ (0, 1].For a given domain Ω ⊂ R d , we define the eigenvalues λ k of the negative Laplacian with periodic boundary conditions as solutions U k (x) of All eigenvalues of Eq. ( 26) are real and non-negative.Then, the solution of Eq. ( 25) can be calculated by hands of a Fourier approach and separation into a homogeneous and inhomogeneous problem.With coefficients c k ∈ R, the solution reads The coefficients can be determined by using the Fourier series of the initial conditions.Analogous formulations are possible for different boundary conditions as well.The solution in Eq. ( 27) is decaying in time towards the ambient temperature T ∞ .
Remark 6.The combustion-free reaction-diffusion problem excluding advection shows a decaying and levelling solution behaviour.
This levelling effect gives the dissipation of the system, compare Eq. ( 5).The convection heat flux Qconv can be integrated over space, resulting in an expression for the energy dissipation rate as The change of E is then given by where we use the Green-Gauß divergence theorem and the zero-flux boundary conditions.The differential equation in E has the solution clearly showing the dissipative nature of the spreading and levelling.The energy dissipation rate tends to zero as the system's temperature approaches the ambient temperature.
Remark 7. The energy dissipation is still given if we replace the constant diffusion by a non-linear diffusion k t (T ) = 4σδT 3 + k.
We explore the energy dissipation now numerically.Fig. 6 shows the solution of Eq. ( 25) for the parameters given in Tab. 1.We impose periodic boundary conditions to ensure that the total energy is not affected by the fluxes across the domain boundaries.Fig. 6(left) shows the decrease of the energy dissipation rate, E, following Eq.( 30) and for comparison and validation of the numerical solver, the value of E computed by Eq. ( 28) for the numerical solution.Simulation snapshots at different time steps are shown in Fig. 6(right).The maximum temperature at every time step is denoted in each of the snapshots.As time evolves, the heat diffuses from an initial maximum temperature in the centre of the domain.The maximum temperature reduces from T = 400 K to almost ambient temperature T = 300.8K after t = 1000 s.Energy dissipation is maximum in the beginning and tends to zero when the solution evolves in time, due to the levelling effect.The energy dissipation rate for the numerical solution fits the analytical estimation, which is a validation for the numerical solver.

Combustion-free advection-diffusion-reaction sub-model
After regarding the energy dissipation of the advection-free model, we focus on the wind-driven propagation of wildfire.In this sub-model, in addition to the energy dissipation and diffusion mechanisms discussed in Sec.3.2.1,temperature is also advected.Combustion is still not regarded, so Y ′ = 0. We consider the partial differential equation for the temperature as with periodic boundary conditions.Assuming a spatially constant wind velocity v, the equation reads and the solution can be calculated according to Sec. 3.2.1 by substituting the space variable with x − v.The energy dissipation in Eq. ( 30) remains identical because the wind only shifts the whole solution in space and the periodic boundary conditions ensure that the total energy is not affected by the fluxes across the domain boundaries.Fig. 7(left) shows the evolution of total energy inside the domain based on the combustion-free ADR sub-model in Eq. ( 31) with periodic boundary conditions and a wind velocity that is pointing towards the direction (1, 1) T .Comparison of the analytically obtained energy dissipation rate (line) and the numerical solution (circles) shows that the numerical solution preserves the energy as expected.Fig. 7(right) shows simulation snapshots at different time steps.The wind advects the high temperature from the centre of the domain at t = 0 to the upper right-hand corner at t = 210 s towards the lower left-hand corner at t = 756 s due to the periodic boundary conditions.Along the way, the temperature is subject to diffusion such that the maximum temperature decreases as the simulation evolves.As in Fig. 6, the maximum temperature decreases until it nearly reaches the ambient temperature after t = 1000 s.
The combustion-free ADR model in Eq. ( 31) does not have travelling wave solutions because the reaction function has only one stationary state for T = T ∞ .The energy is not conserved which is shown by the non-zero energy dissipation rate in Fig. 7(left).

Advection-free sub-models
As a next sub-model, we investigate the model in Eq. ( 14) without wind by setting v = 0.This results in a reaction-diffusion equation coupled to an ordinary differential equation as We impose periodic boundary conditions and set the diffusion to a constant value k t (T ) = k.
The solution of Eq. ( 33) is either dominated by diffusion or has the form of travelling waves.Both solution types are illustrated in Fig. 8 that shows them for some time steps.The diffusion dominated solution (green lines) occurs for parameter choices with comparable small values of H and h that diminish the effect of the reaction terms.The temperature profile diffuses as time progresses.The travelling wave solution (red lines) shows a distinct temperature profile with a steep gradient emerging and travelling through the domain as the simulation evolves.

Travelling wave speed
Analytical results proving upper and lower bounds for the travelling wave speed in the full ADR wildfire model are still an open issue.Here, we investigate the emergence and propagation of travelling waves through numerical simulations and give analytical estimates for the propagation speed.We choose constant initial conditions on a bounded subdomain with the parameters from Tab. 1 and then study variations of the parameters h and k in Tab. 2 and Tab. 3. In particular, we consider the computational domain [0,500] m and we set as initial condition Fig. 9 shows travelling wave solutions that result from the simulation of case A-6 in Tab.  is constant and equal for the spread in both directions.The wave profile of the temperature is steep on the propagation front and less steep where the temperature decays.The biomass burns down to a small amount of remaining biomass and the value is different in the support of the initial condition.Our results are consistent with the findings reported in the literature, for example, [23].
The wave speed and the profile of the wave front depend on the model parameters.Fig. 10(left) shows the influence of the parameter h on the position of the wave front.In all cases, the dependency of the wave front position over time is linear, after a short transient time.The slope of these lines is the wave speed and decreases for increasing values of h, due to a higher heat release into the atmosphere.As seen in Fig. 10, the maximum temperature decreases and the front profile narrows for increasing values of h, as the energy dissipation (i.e.heat release into the atmosphere) enhances.
The observed larger maximum temperature and larger wave speed are connected.For a larger maximum temperature, the spatial gradient of the temperature that drives the wave propagation is larger and thus, leads to a larger wave speed.This observation is as well supported by numerical studies of the influence of the diffusion parameter k in Fig. 11, where the six different sub-cases in Tab. 3 are computed.Here, larger diffusion parameters lead to a higher wave speed and a wider wave profile.Recall that the diffusion term models the heat transfer due to conduction and radiation.
In the literature, results for travelling waves have been obtained for sub-models of the ADR wildfire model.For example, the evolution of a travelling wave solution was shown in [4] for a model that neglects the effects of Newton's cooling law in Eq. ( 5).Further neglecting advection, the speed of a travelling wave is bounded by a sub-solution and a super-solution with c ⋆ < c < c ⋆ .The values for c ⋆ and c ⋆ are given in transformed non-dimensional parameters in [4].A transformation to the dimensional variables and parameters of Eq. ( 14) gives lower and upper bounds as for the travelling wave speed c.The extension of the proof to the full model including the heat exchange with the environment is still an open problem.In the following, we use the values of c ⋆ and c ⋆ for comparison with numerically obtained wave speeds.We apply a different approach by using the system again to a reaction-diffusion equation.In [16], for   reaction-diffusion equations of the form with a monostable reaction function f , the existence of a travelling wave front is proven.The influence of advection on the front speed in case of a monostable reaction-diffusion equation is discussed in [1].
The extension of the results to an advection-diffusion-reaction model coupled with an ordinary differential equation still needs to be addressed and requires similar results for super-and sub-solutions like the proof of the existence of a travelling wave solution.
Therefore, we regard the reaction-diffusion equation for a fixed Y .Then, the conditions of [16] for a travelling wave solution are fulfilled.The reaction function has two roots, where one is at T = T ∞ and the other one is the solution of for fixed Y .Varying Y will vary the root, shifting the value for smaller Y to smaller T .Further, the reaction function f (compare Eq. ( 36)) of Eq. ( 37) is positive between the two roots and the derivative of f is positive for the first and negative for the second root.Then, a travelling wave solution of the single reaction-diffusion equation in Eq. (37) exists and the linearised wave speed of the travelling wave is where f is the reaction-function according to Eq. ( 36).This linearised wave speed may give a lower limit, see [16].Here, it is only derived for a sub-model and we investigate its approximation quality numerically.Fig. 12 compares the bounds in Eq. ( 35) for the model with h = 0 and numerical simulations.Additionally, the linearised minimal wave speed of Eq. ( 38) for the model in Eq. ( 37) is shown.The upper bound for h = 0 in Eq. ( 35) is an upper bound for any h.The dependency of the linearised propagation speed in Eq. ( 38) on h is non-linear with a decrease for increasing h.This behaviour is expected from the results plotted Fig. 10, where we found that the front speed decreases with increasing h.From an application point of view, this insight is useful because even if h cannot be determined precisely, a lower estimate for h gives conservative estimates for the speed of the fire front.However, the lower bound of Eq. ( 35) is not a lower bound of the model with h > 0. For the case h = 0, we find a dependency proportional to √ k of the upper and lower bound of the wave speed on the diffusion parameter k.The upper bound for the wave speed in Eq. ( 35) is again an upper bound for the numeric calculated wave speed.The linearised minimal wave speed from Eq. ( 38) gives a rough estimate for the dependency on k but not exact bounds.
In Fig. 13, we compare the travelling wave speed for different initial biomass values.The smaller the initial biomass is, the slower propagates the combustion wave.This is qualitatively as well the interpretation of the linearised wave speed in Eq. ( 38).The dashed line in Fig. 13, bottom right, gives the dependency of the linearised travelling wave speed depending on Y 0 .The travelling wave speed was gained from fixing the biomass Y and regarding only the reaction-diffusion equation in Eq. (37).This rough approximation of the coupled dynamics shows good results in the computations.
The simulations in Sec.3.3.1 use a constant diffusion coefficient of k t (T ) = k ∈ R. In [7], the influence of different diffusion coefficients k t (T ) was studied numerically.A qualitative comparison of constant diffusion and non-constant diffusion from [7] shows that the burnt area is smaller if a constant diffusion k t (T ) = k is used.

Terminal biomass
The analysis of the lumped sub-model in Eq. ( 16) in Sec.3.1 provides an upper bound for the remaining terminal biomass after the combustion process by the tipping line.A limitation of this analysis is that spatially explicit diffusion processes can not be modelled through the lumped model.Thus, the spatially explicit diffusion mechanism in Eq. ( 33) may lead to a different tipping line and consequently a shift in the upper bound.In spatially explicit sub-models, at any point p in the domain, the diffusion process leads to a rise in temperature before the combustion process reaches it.During the combustion process, heat is transferred to cooler regions by the diffusion mechanism.The increase of temperature at p is amplified by diffusion from warmer regions, leading to a shift of the tipping line in the phase plot to the right-hand side.The diffusion-induced decrease of temperature leads to a shift of the tipping line to the left-hand side in the phase plot.This is illustrated in Fig. 14 that compares trajectories in the phase space of the lumped sub-model in Eq. ( 16) and the advection-free sub-model in Eq. ( 33) for case A-6 in Tab. 2 with the initial condition in Eq. (34).For the latter, trajectories for two different points (p = x) are plotted.One of these points is in the support of the initial condition, where the diffusion effect is small; the second point x further away from the support of the initial condition, where the travelling wave behaviour dominates.Due to the initial temperature distribution, there is a difference in the initial value in the phase space portrait in Fig. 14 for the two points p in space.The trajectory for x = 250 m starts with an initial temperature above the activation temperature, the trajectory for x = 350 m has an initial temperature smaller than the activation temperature.The trajectory for a point where diffusion is small (x = 250 m) follows the trajectories of the ordinary differential equations.In contrast, the trajectory for a point where the travelling wave front passes (x = 350 m) has a decreased temperature compared to the trajectories of the ordinary differential equations.In this particular case, the remaining biomass in both sub-models is similar.However, this is not true in general.In some cases, we observed a decreased biomass for the sub-model in Eq. ( 33) compared to the lumped model in Eq. ( 16).
The influence of the initial biomass on the terminal biomass is qualitatively similar to the lumped model in Fig. 4. Fig. 13 compares the trajectories for different initial values Y 0 .We discussed in Sec.3.3.1 the slowing down of the travelling wave speed by smaller initial biomass.In contrast, the initial biomass has a low impact on the terminal biomass.

Summary of the analysis of the hierarchical sub-models
The spatially lumped sub-model in Sec.3.1 has a tipping line for the shift between heating and cooling which gives an upper limit for the remaining biomass after the combustion process.Analysis of the system reveals a continuous line of stationary states due to the switching function in the reaction function modelling the activation and deactivation of the combustion process.The terminal biomass of the spatially lumped submodel reacts less sensitive to changes in the initial biomass but more sensitive to changes of the activation temperature T ac and changes of the bulk density ρ 0 or the heating value H. Changes in the activation temperature can be interpreted as the influence of management strategies such as pre-emptive watering while changes in the vegetation result in a different bulk density and heating value, in the context of the ADR model herein considered.
We define an energy dissipation functional for the combustion-free sub-model in Sec.3.2 which allows to validate the numerical model for the diffusion and ambient cooling.The energy dissipation does not change if advection caused by wind is introduced.Therefore, the validation of the numerical model is valid for the combustion-free sub-model as well with advection.Meanwhile, this sub-model shows the diffusive behaviour and the exchange with the environment, both leading to a decaying temperature.
The advection-free model in Sec.3.3 including combustion has travelling wave solutions with steep wave profiles.Such travelling wave behaviour is not caused by wind and advection, but it is a result of the interplay between reaction and diffusion.After a short transient phase, a wave profile forms and propagates with a constant speed in space.The wave speed depends on the reaction parameters and the diffusion constant, compare Fig. 10 and 11.As a reference value for the wave speed, we find upper and lower bounds for h = 0.For a constant biomass, the linearised wave speed gives an approximation for the dependency on the wave speed as well.When the initial biomass is reduced, there is a noticeable decrease in wave speed, compare   4 Insights from the full model

One-dimensional wildfire simulation
The full ADR wildfire model allows to study the interplay of combustion, temperature exchange with the environment, diffusion, and advection.Here, we run several one-dimensional simulations using the full ADR wildfire model in Eq. ( 14) with different wind speed v as summarised in Tab. 4. In general, the advection speed of the wildfire should always be lower than the wind speed.In fact, a correction factor must be used to obtain the advection speed [15], but herein, for the sake of simplicity, we will set both of these speeds equal.
Adding wind into the system dynamics changes the amplitude and propagation speed of the waves propagating in the same and in the opposite direction of the wind.Fig. 15 shows the wave profile and the spread of the waves for different times t and a constant wind speed v = 0.08 m/s.Compared to the simulations without wind in Fig. 9, the temperature wave propagating to the left and to the right have different amplitudes, profiles and speed.The biomass consumption differs for the two directions.While the temperature for the wave propagating to the right is much higher, the remaining biomass is higher in this case.This suggests that temperature is not the primary control of the remaining biomass.This observation is consistent with our study of the lumped sub-model in Fig. 4. A higher temperature leads to a shift to the right in the phase space of Fig. 4. The lumped sub-model does not react sensitive to this change with respect to the remaining biomass.Hence, the full model highlights an additional dependency of the remaining biomass.As the wave propagating to the left has a much lower wave speed, the combustion time at a certain position x is much longer than for the wave propagating to the right.This leads to a longer burning process at this position Here, the phase plot shows an increased temperature for the full ADR wildfire model trajectory, compared to the lumped sub-model, which is caused by the transport of temperature in the wind direction.The opposite effect is seen for the green curve, where the wind direction and the wave propagation are opposite.As a consequence, the temperature is lower for medium temperatures but the combustion process at one point in space takes longer.This leads to a lower temperature after the combustion process but to a lower terminal biomass at the equilibrium state.Fig. 16(right) shows the profiles of the waves propagating to the left and to the right.The profile is steeper on the activation front and the amplitude of the temperature is higher for the wave where wind and propagation have the same direction.
Fig. 17 shows the wave profiles of the temperature for different wind speed.The amplitude and the speed of the wave propagating in the opposite direction to the wind decrease with higher wind speed.For a certain wind speed, the opposite wave even stops propagating and only one wave is travelling in the wind direction.This qualitative change of the system's behaviour is as well depicted in Fig. 18.
The propagation speed of the travelling wave in the opposite direction to the wind decreases approximately linearly until a certain threshold is reached.Then, the speed drops to zero-in our case at v = 0.1 m/s.For this speed, which depends on the fixed parameters in Tab. 1, the solution switches from two travelling waves to only one.The travelling wave propagating in the opposite direction of the wind stops.
Fig. 19 compares the evolution of temperature and biomass over time for a wind speed with two travelling waves and with only one.For v = 0.09 m/s, the travelling wave propagating opposite to the wind direction still exists, even if the propagation speed is slow.In contrast, for v = 0.093 m/s, the wave propagating against the wind stops spreading after some time.The temperature profile does not form towards a travelling wave profile and there is only one wave in the direction of the wind, moving to the right-hand side.

Two-dimensional wildfire simulation
We simulate wildfire in a forest ecosystem with a domain of [0, 500]×[0, 500] m 2 .The forest features spatially heterogeneous initial biomass organised in patches, defined as a multi-scale random distribution.Random maps of different scales (i.e. 2, 5, 10 and 100 m) are combined to produce the initial condition for the biomass, Y (x, y, 0).The initial biomass distribution resembles areas with reforestation, for example in areas   Simulation results are plotted in Fig. 20, where at the top, temperature contours and at the bottom, the burned biomass contours at the end of the simulation are overlaid with the terminal and the initial biomass distribution, respectively.The burnt area near to the ignition point preserves its symmetry more or less, but once the fire reaches a patch with low biomass, the symmetry is broken.The wildfire propagation is slowed down in this region.As the wildfire propagates through the domain, encountering more heterogeneity, the symmetry breaks down even further.The expansion towards the north is severely limited by patches of low biomass that draw a "border" along the diagonal of the domain.Besides, when the fire front reaches the patches with lower initial biomass, the temperature of the fire front decreases and thus the severity of the fire, as described in Sec.3.3.2.This is how firebreaks and prescribed burn work in wildfire management.This slower propagation speed in areas with lower biomass was discussed in Sec.3.3.1 for the advection-free model.The linearised wave speed in Eq. (38) depends on the biomass.Fig. 13 also shows this dependency of the advection speed on the biomass.We further observe that the wind speed is sufficiently high to stop the wildfire propagation in the opposite direction.This is consistent with the one-dimensional wildfire model results in Fig. 18.

Conclusions
The advection-diffusion-reaction wildfire model is a non-linear system of several coupled mechanisms between the energy of a fire, represented by its temperature in the context of this model, and the biomass.Sub-models including selected mechanisms form a hierarchical model family and shed light on the complex mechanisms one by one.Numerical simulations of the sub-models support the analytical findings and extend the insights by providing non-formal limit case observations.
In our discussion, we focus on three interlinked properties that we consider important for wildfire management applications: the terminal biomass, the maximum temperature and the propagation speed of the wildfire.The terminal biomass indicates the fraction of the forest ecosystem that went unharmed by the wildfire and is relevant to assess the damage and the potential of recovery, see [13,18].The maximum temperature is an indicator of damage, where higher temperatures will destroy more biomass.The propagation speed of the wildfire influences the amount of burnt area per time.
The spatially lumped ordinary differential equation model in Sec.3.1 provides insight into the relevance of the reaction parameters on the terminal biomass and the maximum temperature.In Sec.3.2 and 3.3, diffusion and advection mechanisms are added to obtain spatially explicit sub-models.
A sensitivity analysis of the lumped ordinary differential equation model shows a high sensitivity of the tipping line for the maximum temperature to the activation temperature, the bulk density and the heating value.This tipping line connects the two relevant aspects of terminal biomass and maximum temperature: shifting the tipping line further up means having a lower maximum temperature and a possibly higher remaining biomass after the fire-see Fig. 5.Meanwhile, the initial biomass has a rather low influence on the terminal biomass-see Fig. 3, Fig. 13 and Fig. 20 for the reaction-diffusion model.In other words, according to the model, reducing the biomass before the occurrence of a wildfire will not help the forest ecosystem to recover more quickly afterwards.However, the analysis of the reaction-diffusion model reveals one advantage of reducing the initial biomass: a lower initial biomass might lead to a reduced wildfire propagation speed-see Fig. 12 and Fig. 13.We found a surprisingly good approximation of the propagation speed by linearisation techniques combined with fixing the biomass Y , compare Eq. (38) and Fig. 13.Diffusion drives the temperature to propagate as a travelling wave with a fixed wave profile.The propagation speed depends on the reaction parameters, the wind speed, and the diffusion parameter.A strong wind affects both the terminal biomass and the wildfire propagation speed.In the direction opposite to the wind, the terminal biomass is smaller and the propagation speed is slower, which might form a natural fire break.In the wind direction, we have a faster wildfire propagation but a higher terminal biomass.The influence of the wind is therefore ambiguous: on the one hand, the wind enforces the wildfire in the wind direction and leads to a much faster propagation and larger burnt area in a shorter time; however, it also has a potential benefit of reducing damage to vegetation, resulting in higher terminal biomass.On the other hand, the area opposite to the wind direction remains unaffected, which serves as a small advantage compared to scenarios with no wind.We observed this model property as well in the two-dimensional simulation in Fig. 20.
The wind-driven advection mechanism in the present ADR wildfire model is a rather simplified representation of real-world processes.In our opinion, it also presents the most interesting area of future research activities.From a mathematical perspective, the shift of the system's behaviour depending on the wind speed is an emergent property that we observe numerically.A full proof and analysis of travelling waves for the full system is still an open problem.Further investigations into the connection between the approximated models and the full model may concretise the quality of the linearised propagation speed.
With regard to model development, the influence of the wind on the fire dynamics could be simulated with higher physical fidelity by coupling the ADR wildfire model to an atmospheric model.This coupling is computationally challenging but promises a more precise prediction of the propagation direction and speed.Such a coupled model in combination with real-world forest structures in a two-dimensional setting would provide the option to test wildfire management strategies like firebreaks or planting formations.It is also crucial for understanding plume-dominated fires, which are expected to become more frequent due to climate change.For the simulation, we choose the following parameters: u = 5, ρ 0 = 1, c = 1, α = k t = 10, R = h = 0.01, C = 0.001 and x 0 = 250.The spatial domain is Ω = [0, 1000] m and the final time is set to t = 100.The solution is computed using N = 100 computational cells and a Courant-Friedrichs-Lewy number of CFL = 0.1.Periodic boundary conditions are set.The numerical solution is computed using a 1-st, 3-rd, 5-th and 7-th order advection scheme.Fig. 23 shows the computed temperature profile at t = 100 s.The lower order schemes (1-st and 3-rd) introduce a high numerical diffusion, which would overestimate the conduction and radiation heat transfer processes.On the other hand, the 5-th and 7-th order schemes provide accurate results for this grid.The results evidence that the advective terms must be computed with high order of accuracy to avoid unphysical diffusion.

A.3 Dependency of the travelling wave speed on the discretisation
The simulation of the reaction-diffusion equation in Eq. (25) shows travelling wave solutions.The wave speed there depends on the number of discretisation points used in the spatial discretisation, see Fig. 24.
In the simulations the parameter h is given by h = 1.The wave speed increases with a larger number of spatial discretisation points with a tendency towards a threshold wave speed.This phenomenon is known for the discretisation of wave-like solutions, compare [35, p. 169f.].

Figure 1 :
Figure 1: Representation of the two-dimensional spatial domain Ω and heat fluxes relevant to wildfire spreading: the conductive flux Qcond , the convection heat flux Qconv , radiation heat flux Qrad , and the reaction heat flux Qreac .The convective flux is driven by the wind velocity v.All heat fluxes can be defined for an arbitrary control volume V ⊂ Ω. Colours inside the domain are meant to represent the temperature distribution.

Theorem 3 .
The solutions of (16) are bounded with (T, Y ) ∈ D = [T low , T max ] × [0, 1] for any initial values (T 0 , Y 0 ) ∈ D and T max sufficiently large.Proof.We start with showing the boundedness of Y .

Figure 2 :
Figure 2: Phase space portrait of the system, showing the heating (red) and cooling (blue) regions of the trajectories, separated by the dashed gray line.Trajectories of the solutions for the initial conditions T 0 = 470 K and Y 0 = {0.2,0.4, 0.6, 0.8, 1.0} are plotted with black continuous lines.

Figure 4 :
Figure 4: Terminal biomass Y ⋆ as a function of the initial temperature and biomass.

Figure 5 :
Figure 5: Sensitivity analysis of the tipping line in Eq. (20).Left: Change of the parameter m = ϱ 0 HA.Right: Change of the activation temperature T ac .

Figure 6 :
Figure 6: Numerical example with the parameters in Tab. 1 for the combustion-free reaction-diffusion model in Eq. (25) with periodic boundary conditions.Left: Evolution of the total energy dissipation rate in the system calculated analytically (line) and numerically (circles).Right: Simulation snapshots of the temperature distribution in the domain at different time instants.Red colour indicates high and yellow colour indicates low temperature.

Figure 7 :
Figure 7: Numerical example with the parameters in Tab. 1 for the combustion-free advection-diffusionreaction model in Eq. (31) with wind velocity in direction (1, 1) T and periodic boundary conditions.Left: Evolution of the total energy dissipation rate in the system calculated analytically (line) and numerically (circles).Right: Simulation snapshots of the temperature distribution in the domain at different time instants.Red colour indicates high and yellow colour indicates low temperature.

Figure 8 :Table 2 :
Figure 8: Comparison of numerically obtained solutions dominated by diffusion (green) and travelling waves (red).Each line represents the solution at a different time instants.The initial condition is represented as black dashed line.The travelling wave solution has been computed using the default parameters in Tab. 1, whereas the diffusion-dominated solution has been computed setting H = 100 kJ/kg and h = 0.01 kWm −3 K −1 .The solutions are plotted every 50 s.
Fig.9shows travelling wave solutions that result from the simulation of case A-6 in Tab. 2 and the parameters from Tab. 1.The initial condition in Eq. (34) evolves to a fire front propagating to the righthand side and a fire front propagating to the left-hand side.The amplitude of the travelling wave solution

Figure 9 :
Figure 9: Numerical example with travelling wave solution for case A-6 in Tab. 2 on a domain x ∈ [0, 500] with N = 2000 cells.All other model parameters are from Tab. 1.The solution is plotted for every 150 s.

Figure 10 :
Figure 10: Numerical example with travelling wave solution to illustrate the sensitivity of the parameter h.The values of h are chosen according to Tab. 2. Left: Position of the wave front depending on the parameter h.Right: Profile of the wave front depending on the parameter h.The parameter k = 2 kW•m −1 •K −1 is fixed.

Figure 11 :
Figure 11: Numerical example with travelling wave solution to illustrate the sensitivity of the parameter k.The values of k are chosen according to Tab. 3. Left: Position of the wave front depending on the diffusion parameter k.Right: Profile of the wave front depending on the parameter k.The parameter h = 1 kW•m −3 •K −1 is fixed.

Figure 12 :
Figure 12: Left: Wave speed for different values of h and fixed k.Right: wave speed for different values of k and fixed h.The bounds from Eq. (35) are depicted in blue, and the linearised wave speed of Eq. (38) is displayed for various fixed values Y .

Figure 14 :
Figure 14: Left: Phase space portrait of the solution for case A-6 at x = 250 and x = 350 m with k = 2 kW•m −1 K −1 , h = 4 kW•m −3 •K −1 .Right: Temperature T and biomass Y over time.The initial condition for the biomass is Y 0 = 1.0.

Figure 15 :
Figure 15: Solutions of the model in Eq. (14) each 150 s for constant diffusion and case C-5 of Tab. 4 with v = 0.08 m/s.The domain is [0,500] m and the discretisation uses N = 2000 cells.The parameters have the default values of Tab. 1 with k = 2 and h = 4.

Fig. 13 .
Fig.13.However, this reduction in initial biomass does not have a significant impact on the terminal biomass, which remains relatively unchanged.

Figure 16 :
Figure 16: Phase space trajectories for three points for case C-5 in Tab. 4. At x = 250 m, the initial condition is non-zero and the behaviour is dominated by the ordinary differential equation.The gray lines give the trajectories of the ordinary differential equation, compare Fig. 2. At x = 200 m, the left travelling wave passes by and at x = 350 m the right travelling wave passes by.

Figure 17 :
Figure 17: Profiles of the travelling temperature waves to the left and to the right for fixed parameters and different wind speed.For v = 0.1 m/s, there is no wave propagating to the left.

Figure 18 :
Figure 18: Left: Position of the wave profiles for different wind speed.Right: Propagation speed of the traveling wave profiled depending on the wind speed.

Figure 19 :
Figure 19: Temperature and biomass evolution for different wind speed v = 0.09 m/s (top) and v = 0.093 m/s (bottom) and the parameters of Tab. 1.

Figure 20 :
Figure 20: Two-dimensional numerical simulation of wildfire propagation through a heterogeneous domain.Top: Temperature contours (every 100 s) and terminal biomass.Bottom: Initial biomass and difference between terminal and initial biomass.Results are shown for 500 s.

Figure 21 :
Figure 21: Numerical solution for T (left) and Y (right) computed by the implicit Euler solver (top) and the implicit RK2 solver (bottom).The reference solution is depicted with black solid line.

Figure 22 :Table 6 :
Figure 22: Computed phase space portrait by the the implicit Euler solver (left) and the implicit RK2 solver (right).The reference solution is depicted with black solid line.

Figure 23 :
Figure 23: Numerical solution computed by a 1-st , 3-rd, 5-th and 7-th order, compared with the exact solution (black solid line).Initial condition depicted with dashed line.

Figure 24 :
Figure 24: Position of the wave front depending on the number of space discretisation points.

Table 5 :
Numerical errors for T and Y measured using the L ∞ norm in the full domain.