On Wind Effects in a Hyperbolic Advection–Reaction–Diffusion Forest Fire Model: Analytical Solutions, Stability and Bifurcation Analysis

1. Introduction

Nonlinear differential equations constitute the backbone of mathematical descriptions for a wide class of complex phenomena observed in natural systems and in human activities [1,2,3,4,5,6,7]. Mathematical modelling based on such equations provides a versatile framework for exploring, interpreting, and forecasting the dynamical behaviour of processes arising in biology and medicine [8,9,10,11], physics and engineering [12,13,14,15,16,17], population dynamics and ecological applications [18,19], and many other scientific fields.

From an analytical perspective, nonlinear models are frequently investigated using constructive techniques for obtaining exact solutions [20,21,22,23,24,25,26,27,28], as well as methods from the qualitative theory of dynamical systems, including stability analysis, bifurcation theory, and the study of complex and chaotic behaviour [8,12,13,14,15,16,17]. In the present study, we combine a dedicated analytical solution technique with qualitative dynamical analysis in order to elucidate the mechanisms governing wind-assisted forest-fire propagation.

Mathematical modelling of forest fire propagation has a long tradition, driven by the need to understand, predict, and control wildfire spread under varying environmental conditions. Wildfire dynamics result from a complex interplay between combustion kinetics, heat transfer, wind forcing, and terrain effects. Consequently, a wide spectrum of modelling approaches has been developed, ranging from empirical and semi-empirical formulations to models based on ordinary differential equations (ODEs) and partial differential equations (PDEs). Comprehensive reviews highlight the diversity of physical, quasi-physical, and phenomenological models proposed in the literature [29,30,31].

Among PDE-based descriptions, parabolic reaction–diffusion and advection–reaction–diffusion equations have become standard tools for investigating wildfire propagation. In many physically motivated formulations, the governing system consists of coupled nonlinear evolution equations for temperature and fuel mass fraction derived from conservation laws and combustion kinetics. Such parabolic models enable systematic investigation of front propagation, wave speeds, and qualitative front structure [32,33,34,35].

In most cases, reaction–diffusion wildfire models employ Fisher–KPP-type nonlinearities to describe the balance between local ignition and fuel depletion through nonlinear growth and saturation mechanisms. However, due to their parabolic character, these models predict infinite propagation speed—an assumption that becomes questionable when modelling real wildfire fronts advancing with finite transport velocities [36].

To overcome this limitation, hyperbolic reaction–diffusion and reaction–transport models have been introduced. By incorporating relaxation (inertial) effects, these formulations ensure finite propagation speed and account for memory effects in the transport process. Hyperbolic corrections substantially modify both the propagation speed and the internal structure of fronts compared to their parabolic counterparts [36,37,38]. These effects are particularly relevant for forest fires, where ignition delay and finite heat-transfer rates play an essential role.

In realistic wildfire scenarios, advection induced by wind or terrain slope decisively influences the directionality, asymmetry, and stability of propagating fire fronts. Wind-driven fire spread has been investigated using a variety of modelling approaches, including operational fire-spread models, cellular automata, level-set techniques, and coupled atmosphere–fire simulations [39,40,41,42,43]. Within PDE frameworks, advection–reaction–diffusion models exhibit strong anisotropy and complex wave morphologies that cannot be captured by purely diffusive formulations [39,40]. When combined with hyperbolic transport, advection produces a rich interplay between reaction kinetics, finite-speed effects, and directional transport [37,38].

The derivation of exact travelling-wave solutions to wildfire models is essential for clarifying fire-front propagation, as such solutions provide a reduced yet physically meaningful description of the front and allow systematic classification of propagation regimes. Analytical insight into these regimes is often complemented by numerical simulations illustrating the influence of parameters such as advection strength, relaxation time, and reaction nonlinearity on front shape and stability [34,35]. In addition, the realistic derivation of exact travelling-wave solutions to nonlinear evolution equations describing real processes is crucial for obtaining reliable information for prediction and control. Such solutions provide benchmark profiles, clarify parameter dependencies, and support further stability and bifurcation analysis. For instance, the qualitative structure of travelling-wave solutions is closely linked to the dynamical properties of equilibria in the reduced travelling-wave systems associated with the investigated PDE models, including their stability type and transitions between different dynamical regimes. These transitions are governed by Hopf bifurcation thresholds marking the onset of oscillatory instabilities. Such mechanisms are well established in the general theory of front propagation into unstable states [44] and in reaction–transport systems with finite propagation speed [37,38].

In this study, we reconsider the hyperbolic reaction–diffusion model for wildfire propagation introduced in [36]. We reformulate the dynamical variable and additionally incorporate an advective transport term to account for wind-induced effects. The resulting two-dimensional planar formulation of the modified model is presented in Section 2. One of the main goals of the present study is to derive exact travelling-wave solutions to the extended model. For this purpose, we apply the Simple Equations Method (SEsM) [45,46], due to its universal character, as emphasized in [47]. A brief description of one of the implemented SEsM (1,1) algorithm variants is also provided in Section 2. In Section 3, physically relevant families of exact solutions of the extended wildfire model are constructed, in accordance with the main dynamical features of the considered process. The steps of the applied SEsM algorithm are formulated as mathematical statements in order to ensure clarity and structural transparency. The obtained analytical solutions are further supported and illustrated by numerical simulations. In Section 4, we perform stability and bifurcation analysis of the equilibrium states of the reduced travelling-wave model and relate the structure of its travelling-wave solutions to the dynamical type of the corresponding equilibria, both analytically and numerically. Finally, in Section 5, the main results are summarized and conclusions are drawn in the context of wind-assisted forest fire propagation.

2. Formulation of the Problem and the SEsM (1,1) Algorithm

2.1. On the Extended Model

We build upon the following hyperbolic reaction–diffusion wildfire model introduced in [36]:

$$\tau n_{tt}+n_t=D{\nabla }^{2}n+F\left(n\right)+\tau F^{\prime }\left(n\right)n_t.$$

(1)

In contrast to the interpretation adopted in [36], where $n(\mathbf{x},t)$ is treated as the density of burning trees, in the present study we interpret $n(\mathbf{x},t)$ as the fire-affected tree fraction (or effective fire–affected tree density).

In Equation (1), the term $n_{tt}$ introduces inertial effects and accounts for finite fire propagation speed, $n_t$ describes the local temporal rate of change of the fire–affected tree fraction, ${\nabla }^{2}n$ represents diffusive interaction between neighbouring regions, and $F\left(n\right)$ is the nonlinear reaction term modeling local ignition dynamics. The coefficient $\tau$ denotes the relaxation time associated with delayed fire flux, while D is the diffusion coefficient corresponding to random fire spread. In [36], the reaction term is given explicitly by

$$F\left(n\right)=rf=r{n}^{\beta }(1-n),\beta \ge 1,$$

(2)

where r is the reaction rate constant (inverse characteristic reaction time), and the parameter $\beta$ quantifies the ignition sensitivity, i.e., the number of fire-affected trees required to ignite a neighbouring unburned tree.

Starting from Equation (1), we incorporate an advective transport term to account for wind-induced effects. Assuming an incompressible wind field, $\nabla ·\mathbf{v}=0$, we have $\nabla ·\left(\mathbf{v}n\right)=\mathbf{v}·\nabla n$. The extended model then reads

$$\tau n_{tt}+n_t=D{\nabla }^{2}n-\nabla ·\left(\mathbf{v}n\right)+F\left(n\right)+\tau F^{\prime }\left(n\right)n_t.$$

(3)

We consider the two-dimensional case $n=n(x,y,t)$. Following the nondimensionalization procedure in [36], we introduce

$$x^{*}={\displaystyle \frac{x}{L}},y^{*}={\displaystyle \frac{y}{L}},t^{*}=rt,{\tau }^{*}=r\tau ,{\mathbf{v}}^{*}={\displaystyle \frac{\mathbf{v}}{\sqrt{Dr}}}.$$

(4)

For simplicity, we drop the asterisks in the sequel, with all quantities understood to be dimensionless.

We focus on the simplest ignition regime $\beta =1$, corresponding to $f\left(n\right)=n(1-n)$. This choice represents a minimal ignition mechanism, where a single fire-affected tree is sufficient to ignite a neighbouring unburned tree. Physically, this regime is appropriate for relatively dense and homogeneous forests, where local heat transfer and flame contact can readily trigger combustion. In this case, fire spread is reaction–limited rather than ignition–threshold-limited.

The final two-dimensional dimensionless extended model takes the form

$$r\tau n_{tt}+\left[1-r\tau \left(1-2n\right)\right]n_t=n_{xx}+n_{yy}-v_0(c_1{n}_x+c_2{n}_y)+n(1-n),$$

(5)

where $v_0$ denotes the magnitude of the dimensionless wind velocity, and $c_1,c_2$ satisfy

$$c_1^2+c_2^2=1,$$

(6)

so that the wind velocity vector is written as $\mathbf{v}=v_0(c_1,c_2)$. The parameter $v_0$ controls the strength of advection, while the unit vector $\mathbf{c}=(c_1,c_2)$ determines its direction.

In the travelling-wave analysis that follows, both possible orientations of the propagating planar front are considered, allowing for downwind and upwind propagation regimes.

2.2. On the SEsM (1,1) Algorithm

Within the SEsM framework, the notation SEsM ($M,N$) indicates that, in the general setting, M nonlinear differential equations are solved using N simple equations. In the present study, we focus on the particular case SEsM (1,1), which is also referred to as the Modified Method of the Simplest Equation (MMSE) [48,49]. This variant has been successfully applied to a wide range of nonlinear problems (For instance, see [50,51,52,53]).

The SEsM(1,1) can be used for obtaining analytical solutions of NPDEs in the following general form:

$$\Omega \left(u\right(x,t),.........)=0$$

(7)

where the left-hand side of Equation (7) is a relationship containing the function u(x, t) and some of its derivatives. A brief outline of the method is provided below.

1.

Introduction of traveiling–wave transformation and construction of the solution of Equation (7) Introducing an appropriate traveiling–wave transformation, the function $u(x,t)=u(\mathbf{k}.\mathbf{x}\pm ct)=u\left(\xi \right)$ is presented in a finite power series:

$$u\left(\xi \right)=\sum _{i=0}^p{a}_{i}f{\left(\xi \right)}^i,$$

(8)

or (as we use in this study)

$$u\left(\xi \right)=\sum _{i=0}^p{a}_i{\left[V_{b_0,b_1,\dots ,b_q}(\xi ;k,l,q)\right]}^i,$$

(9)

where $a_i$$(i=0,1,\dots ,p)$ are constants to be determined, and $V_{b_0,b_1,\dots ,b_q}(\xi ;k,l,q)$ denotes a generalized special function that can present any known exact solution of the chosen simple equation (the ODE). In detail,its explicit form is strongly determined by of the specific form of the simple equation.

2.

Choice of the simple equation The general form of the simple equation can be written as

$${\left({\displaystyle \frac{d^{k}f}{d{\xi }^k}}\right)}^l=\sum _{j=0}^q{b}_j{f}^j\left(\xi \right)$$

(10)

where k is the order of derivative of $f\left(\xi \right)$, l is the degree of derivatives in the defining ODE, q is the highest degree of the polynomials of $f\left(\xi \right)$ in the defining ODE, and $b_j,(j=0...q)$ are the coefficients in the polynomials of $f\left(\xi \right)$, respectively. The general form of Equation (10) allows one to choose a specific simple equation by fixing k, l, q and $b_j$ to values corresponding to an ODE with a known exact solution that matches the expected wave dynamics of the modeled system.

3.

Balance procedure and extraction of the algebraic system Substitution of Equation (8) (or Equation (9)) and the specific form of Equation (10) into Equation (7) leads to a polynomial in powers of of $f\left(\xi \right)$ (or $V_{b_0,b_1,\dots ,b_q}(\xi ;k,l,q)$). The balance procedure considers all possible combinations of powers of the functions above mentioned to establish relations between the degrees of the selected simple equation and the solution (8). One of these balance relations is chosen depending on two basic principles: (1) to ensure that each of the coefficients in front of $f\left(\xi \right)$ (or $V_{b_0,b_1,\dots ,b_q}(\xi ;k,l,q)$) in the polynomial above mentioned has at least two terms; (2) to ensure the accounting of the main physical effects involved in the modeled system. Setting all coefficients in front of a powers of $f\left(\xi \right)$ (or of the specific $V_{b_0,b_1,\dots ,b_q}(\xi ;k,l,q)$) to zero leads to a nonlinear algebraic system that defines the relationships between the coefficients of Equation (8),Equation (10) and Equation (7).

4.

Derivation of exact solutions. Any nontrivial solution of this algebraic system yields an exact analytical solution of the studied nonlinear partial differential equation of kind (8) (or (9)).

3. Derivation of Exact Solutions of Equation (5) Applying SEsM (1,1)

3.1. Travelling-Wave Reduction and Construction of the Solution of Equation (5)

We introduce a planar travelling-wave transformation of the form

$$n(x,y,t)=n\left(\xi \right),\xi ={\ell }_{1}x+{\ell }_{2}y-ct,$$

(11)

where $\ell =({\ell }_1,{\ell }_2)$ is chosen (without loss of generality) as a unit normal to the planar wavefronts $\ell ·\mathbf{x}=\mathrm{const}$, i.e. ${\ell }_1^2+{\ell }_2^2=1$, and $c\in \mathbb{R}$ denotes the oriented travelling-wave speed. The sign of c determines the direction of propagation relative to the chosen orientation of ℓ; reversing the orientation $\ell \mapsto -\ell$ is equivalent to changing the sign of c.

Since the wind velocity $\mathbf{v}=v_0\mathbf{c}$ introduces a preferred direction in the plane, it is natural to parameterize the front orientation relative to $\mathbf{c}$. In two spatial dimensions, the orthonormal basis $\{\mathbf{c},{\mathbf{c}}_{\perp }\}$, where ${\mathbf{c}}_{\perp }=(-c_2,c_1)$, provides a convenient decomposition of any planar propagation direction.

To allow for both orientations of the travelling front, we introduce $\sigma \in \{+1,-1\}$ and set

$$\ell (\mu ,\sigma )=\sigma {\displaystyle \frac{\mathbf{c}+\mu {\mathbf{c}}_{\perp }}{\sqrt{1+{\mu }^2}}},\mu \in \mathbb{R}.$$

(12)

Equivalently,

$${\ell }_1=\sigma {\displaystyle \frac{c_1-\mu c_2}{\sqrt{1+{\mu }^2}}},{\ell }_2=\sigma {\displaystyle \frac{c_2+\mu c_1}{\sqrt{1+{\mu }^2}}}.$$

(13)

Hence,

$${\ell }_1^2+{\ell }_2^2=1,\mathbf{c}·\ell =\sigma {\displaystyle \frac{1}{\sqrt{1+{\mu }^2}}}.$$

(14)

Substituting (11) into (5) and using (14), we obtain the reduced travelling-wave equation

$$\left(1-r\tau c^2\right)n^{\prime \prime }+\left[c-\sigma {\displaystyle \frac{v_0}{\sqrt{1+{\mu }^2}}}-r\tau c\left(1-2n\right)\right]n^{\prime }+n(1-n)=0,$$

(15)

where primes denote differentiation with respect to $\xi$. Only the projection of the wind velocity onto the front normal enters (15), which motivates the definition of the effective (oriented) advection speed

$$v_{\mathrm{eff}}=\mathbf{v}·\ell =\sigma {\displaystyle \frac{v_0}{\sqrt{1+{\mu }^2}}}.$$

Here $\sigma =+1$ corresponds to downwind orientation (ℓ aligned with $\mathbf{c}$), whereas $\sigma =-1$ represents upwind orientation.

The general solution of Equation (5) (or (15)) has the form

$$n\left(\xi \right)=\sum _{i=0}^p{a}_{i}f{\left(\xi \right)}^i,$$

(16)

or (as we use in this study)

$$n\left(\xi \right)=\sum _{i=0}^p{a}_i{\left[V_{b_0,b_1,\dots ,b_q}(\xi ;k,l,q)\right]}^i,$$

(17)

where the value of p is fixed below.

We note that the reaction term $n(1-n)$ possesses two homogeneous equilibria, $n=0$ and $n=1$. In the context of combustion modelling, physically relevant travelling-wave solutions are therefore expected to represent heteroclinic fronts connecting these states, i.e.

$$\underset{\xi \to -\infty }{lim}n\left(\xi \right)=n_{-},\underset{\xi \to +\infty }{lim}n\left(\xi \right)=n_{+},\{n_{-},n_{+}\}=\{0,1\},$$

together with $n^{\prime }\left(\xi \right)\to 0$ as $\xi \to \pm \infty$. This observation guides the subsequent choice of the simple equation within the SEsM framework.

3.2. The Balance Equation and the Algebraic System

3.2.1. The Balance Equation

Lemma 1. 

Let $f\left(\xi \right)$ is the solution of the simplest equation (10) (with $b_q\ne 0$ and $q>1$) in the general solution (16) of (5). Then a consistent dominant–balance analysis of (15), accounting simultaneously for the combined diffusion–relaxation contribution, the wind–driven transport contribution, and the specific nonlinear reaction structure, yields the unique admissible (model-consistent) balance equation leading to a nontrivial polynomial form (16) of $n\left(\xi \right)$ with $p=q-1$.

Proof. 

We substitute (16), together with (10) into (15) and therefore we obtain the dominant powers

$$n\sim f^p,n^{\prime }\sim f^{p+q-1},n^{\prime \prime }\sim f^{p+2q-2},n^2\sim f^{2p},n{n}^{\prime }\sim f^{2p+q-1}.$$

We now perform a dominant-balance comparison term by term in (15), taking into account its physical contributions in (5). Then:

• Combined diffusion–relaxation term versus the reaction term: Balancing $(1-r\tau c^2)n^{\prime \prime }\mathrm{against}n(1-n)=n-n^2$ means balancing the highest power from $n^{\prime \prime }$ with the highest power from $n(1-n)$, i.e. with $n^2$ leading to

  $$p+2q-2=2p⟹p=2q-2.$$
• Combined diffusion–relaxation term versus the linear advection transport part: Balancing $(1-r\tau c^2)n^{\prime \prime }\mathrm{against}\left(c-v_{\mathbf{eff}}-r\tau c\right)n^{\prime }$ yields

  $$p+2q-2=p+q-1⟹q=1.$$
• Combined diffusion–relaxation term versus the relaxation-induced nonlinear transport part: The nonlinear transport contribution is contained in $\left[c-v_{\mathbf{eff}}-r\tau c(1-2n)\right]n^{\prime }=\left(c-v_{\mathbf{eff}}-r\tau c\right)n^{\prime }+2r\tau cn{n}^{\prime }$. Balancing $(1-r\tau c^2)n^{\prime \prime }\mathrm{against}2r\tau cn{n}^{\prime }$ gives

  $$p+2q-2=2p+q-1⟹p=q-1.$$
• Reaction term versus linear advection transport part. Balancing $n(1-n)\mathrm{against}\left(c-v_{\mathbf{eff}}-r\tau c\right)n^{\prime }$ gives (using $n^2$ as the dominant part of $n(1-n)$)

  $$2p=p+q-1⟹p=q-1.$$
• Reaction term versus relaxation-induced nonlinear transport part: Balancing $n(1-n)\mathrm{against}2r\tau cn{n}^{\prime }$ yields

  $$2p=2p+q-1⟹q=1,$$

Collecting the above comparisons, Equation (15) admits several formal balance relations, namely $p=2q-2$, $q=1$, and $p=q-1$. The case $q=1$ leads to a linear simplest equation and is excluded in the present consideration. The balance $p=2q-2$ corresponds to a classical reaction–diffusion scaling, in which the relaxation-induced nonlinear transport contribution $n{n}^{\prime }$ does not enter at dominant order. In contrast, the relation $p=q-1$ is the only balance that preserves, at leading order, the combined diffusion–relaxation second-derivative term together with the relaxation-induced nonlinear transport, while remaining compatible with the wind-driven transport and the nonlinear reaction structure $n(1-n)$. Hence, $p=q-1$ is the unique admissible (model-consistent) balancing condition leading to a nontrivial polynomial form (16) is $p=q-1$. □

The accuracy of the above balance equation will be further confirmed in the next paragraph.

3.2.2. The Algebraic System

Let us fix $q=2$. Then Equation (10) reduces to the generalized ODE of Riccatii:

$${\displaystyle \frac{df}{d\xi }}=b_0+b_{1}f+b_2{f}^2$$

(18)

Then, according the balance equation above defined, Equation (16) reduces to

$$n\left(\xi \right)=a_0+a_{1}f,$$

(19)

and Equation (20) reduces to

$$n\left(\xi \right)=a_0+a_1{V}_{b_0,b_1,b_2}(\xi ;1,1,2)$$

(20)

Lemma 2. 

Let us consider the case where the simple equation used for finding exact solutions of (5) is of kind (18). The application of the SEsM (1,1) with the construction (19) (or (20)) reduces (5) to the following system of nonlinear algebraic equations

$$\begin{array}{ccc}\hfill & & 4a_2^2{b}_{3}cr\tau -8a_2{b}_3^2{c}^{2}r\tau +8a_2{b}_3^2=0\hfill \\ & & 6a_1{a}_2{b}_{3}cr\tau -3a_1{b}_3^2{c}^{2}r\tau +3a_1{b}_3^2+4a_2^2{b}_{2}cr\tau -14a_2{b}_2{b}_3{c}^{2}r\tau +14a_2{b}_2{b}_3=0\hfill \\ & & 4a_0{a}_2{b}_{3}cr\tau +2a_1^2{b}_{3}cr\tau +6a_1{a}_2{b}_{2}cr\tau -2a_1{b}_2{b}_3{c}^{2}r\tau +2a_1{b}_2{b}_3\hfill \\ & & +2a_2^2{b}_{1}cr\tau -12a_2{b}_1{b}_3{c}^{2}r\tau +12a_2{b}_1{b}_3-4a_2{b}_{3}cr\tau +4a_2{b}_{3}c-4a_2{b}_3{v}_{\mathrm{eff}}=0\hfill \\ & & 2a_0{a}_1{b}_{3}cr\tau +2a_0{a}_2{b}_{2}cr\tau +2a_1{a}_2{b}_{1}cr\tau -a_1{b}_1{b}_3{c}^{2}r\tau +a_1{b}_1{b}_3+2a_1{b}_{3}cr\tau -2a_1{b}_{3}c\hfill \\ & & +2a_1{b}_3{v}_{\mathrm{eff}}+2a_2^2{b}_{0}cr\tau -4a_2{b}_0{b}_3{c}^{2}r\tau +4a_2{b}_0{b}_3-2a_2{b}_{2}cr\tau +2a_2{b}_{2}c-2a_2{b}_2{v}_{\mathrm{eff}}=0\hfill \\ & & 2a_0{a}_1{b}_{2}cr\tau -a_0{b}_2{b}_3{c}^{2}r\tau +a_0{b}_2{b}_3+2a_0{b}_{3}cr\tau -2a_0{b}_{3}c+2a_0{b}_3{v}_{\mathrm{eff}}+2a_0{a}_2{b}_{1}cr\tau \hfill \\ & & +a_1^2{b}_{1}cr\tau -a_1{b}_0{b}_3{c}^{2}r\tau +a_1{b}_0{b}_3+2a_1{b}_{1}cr\tau -2a_1{b}_{1}c+2a_1{b}_1{v}_{\mathrm{eff}}\hfill \\ & & +2a_2^2{b}_{0}cr\tau -4a_2{b}_0{b}_2{c}^{2}r\tau +4a_2{b}_0{b}_2-2a_2{b}_{1}cr\tau +2a_2{b}_{1}c-2a_2{b}_1{v}_{\mathrm{eff}}=0\hfill \\ & & a_0^2{b}_{1}cr\tau -a_0{b}_0{b}_2{c}^{2}r\tau +a_0{b}_0{b}_2+2a_0{b}_{1}cr\tau -2a_0{b}_{1}c+2a_0{b}_1{v}_{\mathrm{eff}}+2a_0{a}_1{b}_{1}cr\tau +a_1^2{b}_{0}cr\tau \hfill \\ & & -a_1{b}_0{b}_2{c}^{2}r\tau +a_1{b}_0{b}_2+2a_1{b}_{0}cr\tau -2a_1{b}_{0}c+2a_1{b}_0{v}_{\mathrm{eff}}+a_0-a_0^2=0\hfill \\ & & a_0-a_0^2+a_1{b}_0-2a_2{b}_2^3-27b_3^3+8a_2{b}_2^6-729b_3^6=0\hfill \end{array}$$

(21)

Proof. 

Applying the SEsM $(1,1)$ to Eq. (5) means that we substitute (19) together with (18) and (11)–(14) into (5). Next, we set to zero the coefficients of the resulting polynomials in f in (5). This procedure leads to the algebraic system given in (21). Moreover, it is evident that each equation in (21) consists of at least two terms, which confirms that the balance relation $p=q-1$ has been chosen correctly. Furthermore, exact solutions of (5) (or (15)) can be derived only under this balance equation. □

3.3. Derivation of Exact Solutions of Equation (5)

Several possible analytical solutions of Equation (5) can be derived in dependence on the numerical values of coefficients in Equation (18) leading to use of different special funcion’s forms (For instance, see [28]). Moreover, although the algebraic system admits many non-trivial solutions, we choose a solution in which the coefficients of the simple equation (18) are free parameters. In addition, we ensure that the eventual disappearance of these parameters does not induce algebraic singularities or a collapse of the coefficients in Equation (20). As a result, the same algebraic framework is sufficient to represent all families of solutions associated with the different reductions of the Riccati-type ODE, under small corrections (zeroing) of any of its parameters.

Proposition 1. 

The exact solutions of Equation (5) of the form (20) obtained when Equation (18) is used as a simple equation, which are real–valued and satisfy the physically admissible constraint $0<n\left(\xi \right)<1$, are as follows:

• When $b_0\ne 0,b_1\ne 0,b_2\ne 0$ in Equation (18) and $\Delta =\sqrt{b_1^2-4b_0{b}_2}>0$:

  $$n^{\left(1\right)}\left(\xi \right)={\displaystyle \frac{1}{2}}{\displaystyle \frac{\Delta -b_1}{\Delta }}-{\displaystyle \frac{b_2}{\Delta }}{V}_{b_0,b_1,b_2}^{\left(1\right)}(\xi ;1,1,2)$$

  (22)

  where

  $$\begin{array}{c}\hfill V_{b_0,b_1,b_2}^{\left(1\right)}(\xi ;1,1,2)=-{\displaystyle \frac{b_1}{2b_2}}-{\displaystyle \frac{\Delta }{2b_2}}tanh\left({\displaystyle \frac{\Delta (\xi +{\xi }_0)}{2}}\right)\\ \hfill +{\displaystyle \frac{D}{{cosh}^2\left(\frac{\Delta (\xi +{\xi }_0)}{2}\right)\left(E-\frac{2b_{2}D}{\Delta }tanh\left(\frac{\Delta (\xi +{\xi }_0)}{2}\right)\right)}}\end{array}$$

  (23)

  for $D\ne 0$, $\left|E\right|>\left|{\displaystyle \frac{2b_{2}D}{\Delta }}\right|$, ${\displaystyle \frac{|b_2|}{\Delta }}{\displaystyle \frac{\left|D\right|}{\left|E\right|-\left|\frac{2b_{2}D}{\Delta }\right|}}\le \epsilon ,\epsilon \in \left(0,{\displaystyle \frac{1}{2}}\right]$, where $D,E$ and ${\xi }_0$ are constants and $\xi =l_{1}x+l_{2}y-ct$.

  $$n^{\left(2\right)}\left(\xi \right)={\displaystyle \frac{1}{2}}{\displaystyle \frac{\Delta -b_1}{\Delta }}-{\displaystyle \frac{b_2}{\Delta }}{V}_{b_0,b_1,b_2}^{\left(1\right)}(\xi ;1,1,2)$$

  (24)

  where

  $$V_{b_0,b_1,b_2}^{\left(2\right)}(\xi ;1,1,2)=-{\displaystyle \frac{b_1}{2b_2}}-{\displaystyle \frac{\Delta }{2b_2}}tanh\left({\displaystyle \frac{\Delta (\xi +{\xi }_0)}{2}}\right)$$

  (25)

  for $D=0$ in Equation (23), where ${\xi }_0$ is constant and $\xi =l_{1}x+l_{2}y-ct$.
• When $b_0\ne 0,b_1=0,b_2\ne 0$ in Equation (18):

  $$n^{\left(3\right)}\left(\xi \right)={\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{b_2}{\sqrt{-b_2{b}_0}}}{V}_{b_0,0,b_2}(\xi ;1,1,2)$$

  (26)

  where

  $$V_{b_0,0,b_2}(\xi ;1,1,2)=\sqrt{-{\displaystyle \frac{b_0}{b_2}}}tanh\left(\sqrt{-b_0{b}_2}(\xi +{\xi }_0)\right)$$

  (27)

  for $b_0{b}_2<0$, where ${\xi }_0$ is constant and $\xi =l_{1}x+l_{2}y-ct$.
• When $b_0=0,b_1\ne 0,b_2\ne 0$ and $b_2<0$ and $b_1>0$ in Equation (18):

  $$n^{\left(4\right)}\left(\xi \right)=-{\displaystyle \frac{b_2}{b_1}}{V}_{0,b_1,b_2}(\xi ;1,1,2)$$

  (28)

  where

  $$V_{0,b_1,b_2}(\xi ;1,1,2)={\displaystyle \frac{b_{1}exp\left(b_1(\xi +{\xi }_0)\right)}{1-b_{2}exp\left(b_1(\xi +{\xi }_0)\right)}}$$

  (29)

  where ${\xi }_0$ is constant and $\xi =l_{1}x+l_{2}y-ct$.

Proof. 

Let us consider the case where the simple equations used for finding exact solutions of (5) are of kind (18), where $\Delta =\sqrt{b_1^2-4b_0{b}_2}>0$. We apply the SEsM (1,1) to (5) on the basis of substitutions of (20), (18), and (11)–(14) into (5). As a result, we obtain the system of nonlinear algebraic equations presented in (21). A non–trivial solution of this system is as follows:

$$\begin{array}{ccc}\hfill & & a_0={\displaystyle \frac{1}{2}}{\displaystyle \frac{\sqrt{{b_1}^2-4b_2{b}_0}-b_1}{\sqrt{{b_1}^2-4b_2{b}_0}}},a_1=-{\displaystyle \frac{b_2}{\sqrt{{b_1}^2-4b_2{b}_0}}}\hfill \\ & & r={\displaystyle \frac{\sqrt{{b_1}^2-4b_2{b}_0}}{\tau \left(v_{\mathit{eff}}\sqrt{{b_1}^2-4b_2{b}_0}-1\right)v_{\mathit{eff}}}},c={\displaystyle \frac{v_{\mathit{eff}}\sqrt{{b_1}^2-4b_2{b}_0}-1}{\sqrt{{b_1}^2-4b_2{b}_0}}}\hfill \end{array}$$

(30)

• When $b_0\ne 0$, $b_1\ne 0$, $b_2\ne 0$ in Equation (18): The substitution of $a_0$ and $a_1$ from (30) into (20), together with (23) and (25), leads to the solutions (22) and (24) of (5), respectively. The condition $\Delta =\sqrt{b_1^2-4b_0{b}_2}>0$ guarantees the reality of all coefficients and of the hyperbolic functions appearing in (22)–(25). For $D\ne 0$, the function $V_{b_0,b_1,b_2}^{\left(1\right)}(\xi ;1,1,2)$ contains an additional rational term whose denominator remains nonzero for all $\xi$ provided that

  $$\left|E\right|>\left|{\displaystyle \frac{2b_{2}D}{\Delta }}\right|.$$
  Under this condition, $V^{\left(1\right)}(\xi ;1,1,2)$ is real, smooth and bounded. Moreover, the correction term satisfies

  $$|a_{1}R\left(\xi \right)|\le {\displaystyle \frac{|b_2|}{\Delta }}{\displaystyle \frac{\left|D\right|}{\left|E\right|-\left|\frac{2b_{2}D}{\Delta }\right|}}.$$
  Imposing

  $${\displaystyle \frac{|b_2|}{\Delta }}{\displaystyle \frac{\left|D\right|}{\left|E\right|-\left|\frac{2b_{2}D}{\Delta }\right|}}\le \epsilon ,\epsilon \in \left(0,{\displaystyle \frac{1}{2}}\right],$$

  ensures that the solution $n^{\left(1\right)}\left(\xi \right)$ remains strictly within the interval $(0,1)$ for all $\xi$. For $D=0$, the function $V^{\left(1\right)}$ reduces to the bounded hyperbolic solution $V^{\left(2\right)}$, which yields

  $$n^{\left(2\right)}\left(\xi \right)={\displaystyle \frac{1}{2}}\left(1+tanh\left({\displaystyle \frac{\Delta (\xi +{\xi }_0)}{2}}\right)\right).$$
  Since $tanh(·)\in (-1,1)$ for all finite arguments, it follows that $0<n^{\left(2\right)}\left(\xi \right)<1$ for all $\xi$, with the limiting values attained only asymptotically as $\xi \to \pm \infty$. Hence, no additional constraints beyond $\Delta >0$ are required in this case.
• When $b_0\ne 0$, $b_1=0$, $b_2\ne 0$ in Equation (18): Setting $b_1=0$ in Equation (30) reduces it to

  $$a_0={\displaystyle \frac{1}{2}},a_1=-{\displaystyle \frac{1}{2}}{\displaystyle \frac{b_2}{\sqrt{-b_2{b}_0}}},r={\displaystyle \frac{2\sqrt{-b_2{b}_0}}{\tau \left(2v_{\mathit{eff}}\sqrt{-b_2{b}_0}-1\right)v_{\mathit{eff}}}},c={\displaystyle \frac{2v_{\mathit{eff}}\sqrt{-b_2{b}_0}-1}{2\sqrt{-b_2{b}_0}}}.$$

  (31)
  The substitution of $a_0$ and $a_1$ from (31) into (20), together with (27), leads to the solution (26) of (5). The condition $b_0{b}_2<0$ ensures the reality of $\sqrt{-b_0{b}_2}$ and hence of the function $V_{b_0,0,b_2}(\xi ;1,1,2)$ defined in (27). Since $tanh(·)\in (-1,1)$, this function is bounded for all $\xi$, which implies that $n^{\left(3\right)}\left(\xi \right)$ is real–valued and takes values strictly in $(0,1)$ for all finite $\xi$.
• When $b_0=0$, $b_1\ne 0$, $b_2\ne 0$ in Equation (18): Substituting $b_0=0$ into Equation (30) yields

  $$a_0=0,a_1=-{\displaystyle \frac{b_2}{b_1}},r={\displaystyle \frac{b_1}{\tau \left(v_{\mathit{eff}}{b}_1-1\right)v_{\mathit{eff}}}},c={\displaystyle \frac{v_{\mathit{eff}}{b}_1-1}{b_1}}.$$

  (32)
  The substitution of $a_0$ and $a_1$ from (32) into (20), together with (29), leads to the solution (28) of (5). In order for the solution $n^{\left(4\right)}\left(\xi \right)$ to be real–valued and satisfy $0<n\left(\xi \right)<1$, the parameters must satisfy $b_1>0$ and $b_2<0$. Under these conditions, $V_{0,b_1,b_2}(\xi ;1,1,2)>0$ for all $\xi$, while the prefactor $-b_2/b_1$ is positive, ensuring that $n^{\left(4\right)}\left(\xi \right)$ remains strictly between 0 and 1. There is also alternative choice $b_1<0$ and $b_2>0$ in (18), when $b_0=0$, as is shown in [28], but it leads to negative or unbounded values of $n\left(\xi \right)$ for finite $\xi$ and is therefore excluded as physically inadmissible for the considered model.

□

The analytical calculations and the numerical simulations reported in the following sections were carried out using Mathematica 10.1.0 Standard Edition (Wolfram Research, Inc., USA).

3.4. Numerical Comparative Analysis of Travelling Fronts Based on the Exact Solutions

The exact travelling-wave solutions derived above form a parametric family of monotone travelling fronts. Although all solutions satisfy the same asymptotic constraints and arise from the same travelling-wave reduction, their internal structure depends on the specific form of the simple equation as well as on the selected parameter set. To illustrate the physical implications of these differences in the context of wildfire-front propagation, we compare numerically the travelling-wave profiles given by Equations (22)–(29) in Figure 1. The parameter values used in the numerical simulations for $n^{\left(1\right)}\left(\xi \right)$–$n^{\left(4\right)}\left(\xi \right)$ are listed in Table 1. In Figure 1, all profiles are centred by imposing $n\left(\overline{\xi }\right)=0.5$.

Figure 1 highlights the quantitative distinctions between the profiles. While their asymptotic limits coincide, the spatial rate at which the transition occurs varies significantly. In particular, the wave profile $n^{\left(3\right)}\left(\xi \right)$ exhibits the most abrupt change around the inflection point, indicating the thinnest transition layer. The solution $n^{\left(1\right)}\left(\xi \right)$ (with $D\ne 0$) produces a comparably steep front, though with a slightly smoother central region. The case $D=0$ in $n^{\left(2\right)}\left(\xi \right)$ leads to a visibly wider transition zone, whereas the logistic profile $n^{\left(4\right)}\left(\xi \right)$ displays the most gradual spatial variation, characterised by a smaller maximum gradient and slower approach to saturation. These differences demonstrate that, even within the same travelling-wave framework, the analytical construction influences the effective front thickness and the spatial distribution of the transition dynamics.

In the wildfire context, these structural differences correspond to variations in the sharpness of the combustion interface. Steeper fronts represent thinner reaction zones and more abrupt transitions between unburned and fire-affected regions, whereas broader profiles reflect stronger diffusive smoothing. Thus, the presented family of exact solutions, based on different Riccati-ODE solution structures, provides additional flexibility, enabling systematic control of both front steepness and internal morphology without altering the propagation speed. This flexibility is essential for modelling different fire-spread regimes under varying environmental conditions.

In the presence of wind, the effective advection velocity entering the travelling-wave equation (15) depends on both the wind intensity and the front orientation which is incorporated through the effective velocity $v_{\mathrm{eff}}$. By fixing the wind magnitude and varying the orientation parameter, one isolates the effect of the front inclination on the travelling-wave structure.

The travelling-wave profiles shown in Figure 2 correspond to the exact solution $n^{\left(1\right)}$ (see Eq. (22)), evaluated for several values of the effective advection velocity $v_{\mathrm{eff}}$. Although $v_{\mathrm{eff}}$ does not appear explicitly in the expression (22), its influence enters through the discriminant $\Delta =\sqrt{b_1^2-4b_0{b}_2}$. Thereby, using the relation given in Equation (30), we rewrite $\Delta$ in terms of $r\tau$ and $v_{\mathrm{eff}}$, which yields $\Delta \left(v_{\mathrm{eff}}\right)={\displaystyle \frac{r\tau v_{\mathrm{eff}}}{v_{\mathrm{eff}}^{2}r\tau -1}}$. This expression is well defined provided that $v_{\mathrm{eff}}\ne 0,v_{\mathrm{eff}}^{2}r\tau \ne 1$, which ensure the non-singularity of the travelling-wave solution. Under these assumptions, Figure 2 presents the profiles of Equation (22) for physical parameters $r=1$ and $\tau =0.5$, while the parameters of the solution of an ODE of Riccati are $b_0=0.6,b_1=0.8,b_2=-1,D=0.15$, and $E=0.80$. All profiles are centred by imposing the condition $n^{\left(1\right)}\left(\overline{\xi }\right)=0.5$.

As indicated by the values of $\Delta$ reported in the legend of Figure 2, smaller values of $v_{\mathrm{eff}}$ correspond to larger values of $\left|\Delta \right|$, whereas increasing $v_{\mathrm{eff}}$ leads to a reduction in $\left|\Delta \right|$. Since the dominant contribution to the front slope is governed by $\Delta$, larger $\left|\Delta \right|$ produce sharper and more localised transition layers, while smaller $\left|\Delta \right|$ generate broader and smoother profiles. This trend is clearly visible in Figure 2: decreasing $v_{\mathrm{eff}}$ yields steeper and thinner fire fronts, whereas increasing $v_{\mathrm{eff}}$ results in wider transition regions.Owing to the presence of the deformation term ($D\ne 0$), the effect of $v_{\mathrm{eff}}$ is not merely a uniform horizontal rescaling of the profile. Instead, variations in $v_{\mathrm{eff}}$ alter the effective front thickness while preserving the monotone heteroclinic structure of the solution.

From a wildfire-propagation perspective, this demonstrates that wind-induced advection, mediated by front orientation, directly influences the sharpness of the combustion front: weaker effective advection produces concentrated and abrupt fire lines, whereas stronger effective advection leads to more spatially distributed and smoother burning transitions.

Finally, Figure 3 presents three-dimensional visualisations of the solution $n^{\left(1\right)}\left(\xi \right)$, where $\xi =l_{1}x+l_{2}y-ct$, at three distinct time instants. The numerical simulations are performed for the parameter values $b_0=0.2$, $b_1=0.6$, $b_2=-1$, $D=0.15$, and $E=0.7$. The travelling-wave direction is determined by the unit vector $(l_1,l_2)=(0.8,0.6)$, satisfying $l_1^2+l_2^2=1$, and the wave speed is set to $c=0.6$.As it is seen from Figure 3, at the initial time, the solution exhibits the characteristic sigmoidal structure of a travelling front, smoothly connecting two distinct asymptotic states. As time increases, the entire surface translates rigidly along the direction $(l_1,l_2)$, while its shape, amplitude, and steepness remain unchanged. No deformation or spreading of the transition layer is observed. This confirms that the analytical construction indeed represents a coherent planar travelling wave propagating with constant speed c, consistently with the travelling-coordinate reduction.

From a dynamical perspective, the moving interface corresponds to a connection between two equilibrium states of the reduced travelling-wave system. The invariance of the profile in the moving frame indicates that these equilibria govern the asymptotic states ahead of and behind the front. Therefore, the qualitative behaviour observed in the two-dimensional visualizations naturally leads to a local stability analysis of the corresponding equilibrium points. Such an analysis allows one to determine their dynamical character and to classify the possible propagation regimes of the travelling fronts.

4. Qualitative Analysis of the Equilibria of the Traveling–Wave Equation (15)

For our convenience, we introduce the following notation in Equation (15):

$$A=1-r\tau c^2,B\left(n\right)=c-v_{\mathrm{eff}}-r\tau c(1-2n),f\left(n\right)=n(1-n),$$

(33)

so that it takes the compact form

$$A{n}^{\prime \prime }+B\left(n\right)n^{\prime }+f\left(n\right)=0.$$

(34)

Introducing

$$u_1=n,u_2=n^{\prime },$$

(35)

Equation (34) can be rewritten as the planar system

$$\begin{array}{c}\hfill u_1^{\prime }=u_2,\\ \hfill u_2^{\prime }=-{\displaystyle \frac{B\left(u_1\right)u_2+f\left(u_1\right)}{A}},\end{array}$$

(36)

Below we analyse the equilibria and their local stability in the phase plane.

4.1. Equilibria

Lemma 3. 

The system (36) possesses exactly two equilibria,

$$E_0=(0,0),E_1=(1,0).$$

(37)

Proof. 

Equilibria satisfy $u_1^{\prime }=u_2^{\prime }=0$. From $u_1^{\prime }=u_2$ it follows that $u_2^{*}=0$. Substituting into the second equation yields $f\left(u_1^{*}\right)=0$. Since $f\left(n\right)=n(1-n)$, we obtain $u_1^{*}\in \{0,1\}$, which completes the proof. □

Thus, travelling-wave solutions presented in the previous section correspond to phase-plane trajectories of (36) connecting or returning to the equilibria $E_0$ and $E_1$.

4.2. Local Stability Classification of the Equilibria

Proposition 2. 

Let $A=1-r\tau c^2\ne 0$ and define $B\left(0\right)=c-v_{\mathrm{eff}}-r\tau c,B\left(1\right)=c-v_{\mathrm{eff}}+r\tau c$ with $v_{\mathrm{eff}}=\sigma v_0/\sqrt{1+{\mu }^2}$. The equilibria $E_0=(0,0)$ and $E_1=(1,0)$ of the travelling-wave system (36) are qualitatively classified as follows.

(a)

For $E_0$:

• If $A<0$, then $E_0$ is a saddle.
• If $A>0$ and $B\left(0\right)>0$, then $E_0$ is asymptotically stable; moreover, it is a stable node if $B{\left(0\right)}^2\ge 4A$ and a stable focus if $B{\left(0\right)}^2<4A$.
• If $A>0$ and $B\left(0\right)<0$, then $E_0$ is unstable; moreover, it is an unstable node if $B{\left(0\right)}^2\ge 4A$ and an unstable focus if $B{\left(0\right)}^2<4A$.

(b)

For $E_1$:

• If $A>0$, then $E_1$ is a saddle.
• If $A<0$ and $B\left(1\right)<0$, then $E_1$ is asymptotically stable; moreover, it is a stable node if $B{\left(1\right)}^2\ge -4A$ and a stable focus if $B{\left(1\right)}^2<-4A$.
• If $A<0$ and $B\left(1\right)>0$, then $E_1$ is unstable; moreover, it is an unstable node if $B{\left(1\right)}^2\ge -4A$ and an unstable focus if $B{\left(1\right)}^2<-4A$.

Proof. 

Linearising (36) at an equilibrium $(u_1^{*},0)$ yields the Jacobian

$$J(u_1^{*},0)=\left(\begin{array}{cc}0& 1\\ -{\displaystyle \frac{f^{\prime }\left(u_1^{*}\right)}{A}}& -{\displaystyle \frac{B\left(u_1^{*}\right)}{A}}\end{array}\right),f\left(n\right)=n(1-n),$$

with characteristic polynomial

$${\lambda }^2+T_j\lambda +D_j=0,$$

where

$$T_j=-\mathrm{tr}J\left(E_j\right)={\displaystyle \frac{B\left(u_1^{*}\right)}{A}},D_j=detJ\left(E_j\right)={\displaystyle \frac{f^{\prime }\left(u_1^{*}\right)}{A}}$$

are the Routh–Hurwitz coefficients.

Since $f^{\prime }\left(0\right)=1$ and $f^{\prime }\left(1\right)=-1$, we obtain

$$T_0={\displaystyle \frac{B\left(0\right)}{A}},D_0={\displaystyle \frac{1}{A}},T_1={\displaystyle \frac{B\left(1\right)}{A}},D_1=-{\displaystyle \frac{1}{A}}.$$

Hence $D_0<0$ if and only if $A<0$, implying that $E_0$ is a saddle, whereas $D_1<0$ if and only if $A>0$, implying that $E_1$ is a saddle. If $D_j>0$, stability is governed by the Routh–Hurwitz condition $T_j>0$, which yields the stated conditions on $B\left(0\right)$ and $B\left(1\right)$.

Finally, the node–focus transition is determined by the discriminant

$${\Delta }_j=T_j^2-4D_j,$$

that is,

$${\Delta }_0={\displaystyle \frac{B{\left(0\right)}^2-4A}{A^2}},{\Delta }_1={\displaystyle \frac{B{\left(1\right)}^2+4A}{A^2}}.$$

Therefore,

(a)

For $E_0$ (thus $A>0$),

$$\mathrm{node}\iff {\Delta }_0>0\iff B{\left(0\right)}^2>4A,\mathrm{focus}\iff {\Delta }_0<0\iff B{\left(0\right)}^2<4A,$$

with the transition (double real eigenvalue) at

$$B{\left(0\right)}^2=4A\iff v_{\mathrm{eff}}=c(1-r\tau )\pm 2\sqrt{A}.$$

(b)

For $E_1$ (thus $A<0$),

$$\mathrm{node}\iff {\Delta }_1>0\iff B{\left(1\right)}^2>-4A,\mathrm{focus}\iff {\Delta }_1<0\iff B{\left(1\right)}^2<-4A,$$

with the transition at

$$B{\left(1\right)}^2=-4A\iff v_{\mathrm{eff}}=c(1+r\tau )\pm 2\sqrt{-A}.$$

□

The explicit appearance of $v_{\mathrm{eff}}$ in the above relations is included to facilitate the bifurcation analysis presented in the next subsection.

Remark 1. 

Equation (36) provides a phase-plane representation for the spatial profiles of travelling-wave solutions. In particular, the asymptotic behaviour of a travelling-wave profile near an equilibrium $E_j$ is governed by the linear type of $E_j$, as classified in Proposition 2.

If $E_j$ is a node (stable or unstable), the travelling-wave profile approaches the equilibrium in a monotone manner, and the corresponding wave tail is smooth and non-oscillatory.

If $E_j$ is a focus (stable or unstable), the travelling-wave profile approaches the equilibrium in an oscillatory manner, leading to a wave tail with damped (or growing) spatial oscillations.

The node–focus transition described above reflects a change in the linear classification of a hyperbolic equilibrium and represents a qualitative variation within the same dynamical regime. As such, it does not constitute a bifurcation.

In contrast, transitions between genuinely different qualitative regimes of travelling-wave behaviour, such as the emergence of periodic orbits in the phase plane, require a loss of hyperbolicity of the equilibrium and are associated with true bifurcations.These phenomena will be analysed in the next subsection in the context of Hopf bifurcations with respect to the parameter $v_{\mathrm{eff}}$.

4.3. Bifurcation Analysis with Respect to $v_{\mathrm{eff}}$

In this paragraph we investigate how the qualitative behaviour of the travelling–wave system (36) changes as the effective velocity $v_{\mathrm{eff}}$ varies. Throughout, we use the notation and the local stability classification of the equilibria $E_0$ and $E_1$ established in the previous section.

Proposition 3. 

Assume $A\ne 0$. Then

(a)

If $A>0$, then the equilibrium $E_0$ has a Hopf threshold at

$$v_{\mathrm{eff}}^{H,0}=c(1-r\tau ).$$

At $v_{\mathrm{eff}}=v_{\mathrm{eff}}^{H,0}$ one has $T_0=0$, $D_0>0$ and the eigenvalues are ${\lambda }_{1,2}=\pm i{\omega }_0$ with ${\omega }_0=1/\sqrt{A}$. Moreover,

$${\displaystyle \frac{d{T}_0}{d{v}_{\mathrm{eff}}}}=-{\displaystyle \frac{1}{A}}\ne 0,$$

hence the crossing is transversal.

(b)

If $A<0$, then the equilibrium $E_1$ has a Hopf threshold at

$$v_{\mathrm{eff}}^{H,1}=c(1+r\tau ).$$

At $v_{\mathrm{eff}}=v_{\mathrm{eff}}^{H,1}$ one has $T_1=0$, $D_1>0$ and ${\lambda }_{1,2}=\pm i{\omega }_1$ with ${\omega }_1=1/\sqrt{-A}$, and

$${\displaystyle \frac{d{T}_1}{d{v}_{\mathrm{eff}}}}=-{\displaystyle \frac{1}{A}}\ne 0.$$

Proof. 

At $E_j$ a Hopf threshold (purely imaginary pair in the linearisation) occurs when $T_j=0$ and $D_j>0$ in ${\lambda }^2+T_j\lambda +D_j=0$. For $E_0$ we have $D_0>0\iff A>0$, and $T_0=0\iff B\left(0\right)=0\iff v_{\mathrm{eff}}=c(1-r\tau )$. The eigenvalues are $\pm i\sqrt{D_0}=\pm i/\sqrt{A}$. Transversality follows from $d{T}_0/d{v}_{\mathrm{eff}}=-1/A\ne 0$. The proof for $E_1$ is analogous using $D_1>0\iff A<0$ and $T_1=0\iff B\left(1\right)=0$. □

Remark 2. 

Since ${\lambda }_{1,2}=-{\displaystyle \frac{T_j}{2}}\pm {\displaystyle \frac{1}{2}}\sqrt{T_j^2-4D_j}$, for $D_j>0$ the equilibrium is linearly asymptotically stable if $T_j>0$ and unstable if $T_j<0$. Hence, for $A>0$ the equilibrium $E_0$ changes from stable to unstable as $v_{\mathrm{eff}}$ increases through $v_{\mathrm{eff}}^{H,0}$ (because $B\left(0\right)$ decreases with $v_{\mathrm{eff}}$), whereas for $A<0$ the equilibrium $E_1$ changes from unstable to stable as $v_{\mathrm{eff}}$ increases through $v_{\mathrm{eff}}^{H,1}$.

A genuine Hopf bifurcation (i.e., the birth of a small-amplitude periodic orbit) additionally requires the standard non-degeneracy condition, namely a nonzero first Lyapunov coefficient. The explicit computation of the first Lyapunov coefficient, and thus the determination of the criticality (supercritical or subcritical) of the Hopf bifurcation, is beyond the scope of the present study.

The local stability classification and the bifurcation analysis with respect to $v_{\mathbf{eff}}$ can now be synthesised into a unified regime structure, summarised in Table 2, which organises the qualitative travelling-wave behaviour according to the dynamical type of the equilibria.

4.4. Complete Regime Map, Representive Phase-Plane Portraits and Physical Interpretation

Based on the qualitative classification summarised in Table 2, a complete regime map in the $(c,v_{\mathrm{eff}})$–parameter plane is constructed and displayed in Figure 4.

The regime map is performed for fixed parameter values $r=1$ and $\tau =0.5$, thereby restricting the study to a representative cross-section of the full parameter space. The parameters c and $v_{\mathrm{eff}}$ are considered with both positive and negative signs because they represent oriented quantities. The sign of c determines the direction of front propagation relative to the chosen normal vector, while $v_{\mathrm{eff}}=\mathbf{v}·\ell$ is the signed projection of the wind velocity onto the front normal, allowing for both downwind and upwind regimes. In Figure 4, the advection velocity $v_{\mathrm{eff}}$ is restricted in the interval $[-2.5,2.5]$, while the travelling-wave speed c is varied within $[-1.5,1.5]$. These bounds are chosen consistently with the constraints introduced during the construction of the model and its reduction to a travelling-wave formulation, where the admissible parameter values follow from the imposed physical and mathematical assumptions. In particular, the selected interval for c includes the critical threshold defined by $1-r\tau c^2=0$, while the range of $v_{\mathrm{eff}}$ covers dynamically relevant regimes in which advection varies from weak to dominant relative to diffusion and reaction. As is seen, the regime map visualizes the regions in parameter space where different types of the two equilibrium states $E_0$ and $E_1$ (stable and unstable nodes and foci) are realized, as well as the loci of Hopf bifurcations. In this way, it provides a comprehensive picture of the qualitative dynamics of the reduced traveling–wave system (36) for the specified values of r and $\tau$. We note that varying the values of r and $\tau$ would mainly lead to a rescaling of the regime boundaries and a redistribution of the corresponding parameter intervals, while preserving the overall qualitative structure of the dynamical picture.Although the regime map represents the complete classification of all possible dynamical regimes in the $(c,v_{\mathrm{eff}})$ parameter plane, additional restrictions arise from the specific form of the coefficients of the exact travelling–wave solutions of Equation (15) presented in Section 3. In particular, for the selected set of coefficients obtained from the algebraic system (21) (For instance, see (30)–(32)), one derives the relation $c={\displaystyle \frac{1}{r\tau v_{\mathrm{eff}}}}$. Since $r\tau >0$, it follows that $c{v}_{\mathrm{eff}}>0$, that is, c and $v_{\mathrm{eff}}$ must have the same sign. We emphasize that this constraint is not intrinsic to the full dynamical system but is specific to the chosen family of exact solutions corresponding to coefficient sets (30)–(32). For other alternative solutions of the algebraic system (21), this restriction may no hold. Nevertheless, in the representative phase portraits presented below, this analytically induced constraint has been taken into account by selecting parameter pairs $(c,v_{\mathrm{eff}})$ consistent with $c{v}_{\mathrm{eff}}>0$ within the regime map.

To illustrate the dynamical regimes identified in the regime map, we present representative phase-plane portraits together with the corresponding travelling-wave profiles in Figure 5 and Figure 6. Each panel shows trajectories of Equation (36) in the $(n,n^{\prime })$ plane, highlighting the type and stability of the equilibria $E_0=(0,0)$ and $E_1=(1,0)$, as well as the behaviour of the traveling wave near them.

The qualitative phase-plane analysis presented here aims to clarifies under which parameter regimes the analytically constructed travelling-wave solutions correspond to physically meaningful invasion fronts. Since the variable n represents the fire-affected fraction, the equilibria $E_0$ and $E_1$ correspond to the completely fire–unaffected ($n=0$) and fully fire–affected ($n=1$) states, respectively. It is obvious that the travelling-wave front corresponds to a heteroclinic connection in the reduced system (36). In the travelling coordinate $\xi$, such a front can occur only when one equilibrium is a saddle (providing a one-dimensional unstable manifold) and the other is asymptotically stable (providing a stable manifold). Therefore, Figure 5 and Figure 6 illustrate in detail the two physically plausible scenarios for the propagation of the forest fire wave front in accordance with the model system (36).

In Figure 5, the scenario $A=1-r\tau c^2>0$ is illustrated when $E_1$ is a saddle and $E_0$ representing the fire–unaffected (or extinguished) state is a non-saddle. Then the heteroclinic relation $E_1\to E_0$ describes a relaxation towards an unburned forest. If $E_0$ is a stable node (for $v_{\mathrm{eff}}<c(1-r\tau )-2\sqrt{A}$), the disturbances decay monotonically (see Figure 5a). This means that any local ignition will quickly decay and cannot develop into a sustained combustion. This corresponds to a fire-resistant environment, for example due to high moisture content, low ambient temperature, weak reaction intensity or strong heat dissipation. If $E_0$ is a stable focus (for $c(1-r\tau )-2\sqrt{A}<v_{\mathrm{eff}}<c(1-r\tau )$), the disturbances decay oscillatorily (see Figure 5b). Physically, small attempts at ignition may lead to temporary fluctuations in combustion or minor attempts at re-ignition, but the system eventually returns to an unburned state, i.e. the forest remains still stable against ignition. At the Hopf threshold of $E_0$ (for $v_{\mathrm{eff}}=c(1-r\tau )$) the system reaches the ignition limit and prolonged combustion oscillations can occur ((see Figure 5c). After this point, when $E_0$ becomes unstable (for $c(1-r\tau )<v_{\mathrm{eff}}<c(1-r\tau )+2\sqrt{A}$ ), small perturbations grow instead of decaying (see Figure 5d). Physically, this corresponds to ignition instability or self-ignition, i.e. spontaneous sustained combustion can occur. In this regime, the initial heteroclinic relationship with $E_0$ disappears and the system instead evolves to an active combustion state, potentially generating a propagating forest fire front.

In Figure 6 the oposite scenario $A=1-r\tau c^2<0$ is illustrated when $E_0$ is a saddle and $E_1$ corresponding to the fire–affected state, is non-saddle. Then the heterocline is $E_0\to E_1$ and it describes a self-sustained wildfire front. If $E_1$ is a stable node (for $v_{\mathrm{eff}}>c(1+r\tau )+2\sqrt{-A}$ ), the approach along the heteroclinic orbit is monotone (see Figure 6a). The corresponding wave profile exhibits a smooth, non-oscillatory transition from fire–unaffected to fire–affected state. Physically, this represents a steady and robust combustion regime, i.e. the fire propagates with constant structure, the fire–affected fraction stabilizes behind the front, and no secondary fluctuations occur. This regime corresponds to sufficiently dry fuel, strong exothermic reaction, and effective heat transfer, ensuring a stable and continuous wildfire propagation. If $E_1$ is a stable focus (for $c(1+r\tau )<v_{\mathrm{eff}}<c(1+r\tau )+2\sqrt{-A}$ ), the heteroclinic orbit approaches it in a spiral manner (see Figure 6b ). The wave profile then exhibits damped oscillations behind the front. Physically, this reflects transient fluctuations in the fire–affected region due to competition between heat production and dissipation. The fire remains self-sustained, but the post-front dynamics may involve small-scale secondary effects. Although dynamically richer, this regime still corresponds to stable wildfire propagation, since perturbations decay asymptotically. At the Hopf threshold (for $v_{\mathrm{eff}}=c(1+r\tau )$ ), the fire–affected equilibrium loses hyperbolicity and the system reaches marginal stability (see Figure 6c ). This marks the transition to a pulsating combustion regime. Once $E_1$ becomes unstable (unstable focus) (for $c(1+r\tau )-2\sqrt{-A}<v_{\mathrm{eff}}<c(1+r\tau )$ ), its stable manifold disappears and the unstable manifold of the saddle can no longer connect to it (see Figure 6d ). Consequently, the heteroclinic orbit, and thus the traveling-wave front ceases to exist. Physically, this means that the fire–affected state cannot remain stable, i.e. either the fire loses its steady self-sustained character and the front breaks down, or the system transitions to a dynamically unstable, oscillatory wildfire regime with possible front fragmentation and spatio-temporal instabilities.

In summary, from both mathematical and physical perspectives, the existence of a fire travelling-wave front, and equivalently of the analytical travelling-wave solutions derived in Section 3, is directly determined by the stability properties of the equilibria of the reduced system (36). A travelling-wave solution can exist only when the terminal equilibrium ($E_0$ or $E_1$) of the corresponding heteroclinic orbit is asymptotically stable, ensuring the persistence of the associated stable manifold in phase space. Physically, this means that a combustion front can propagate in a sustained manner only if the state behind the front (fire–affected or fire–unaffected) is dynamically stable with respect to small perturbations. The node–focus transition influences only the manner in which the wave tail approaches the equilibrium, either monotonically or through damped spatial oscillations, without affecting the existence of the front itself.

However, once the equilibrium loses stability, typically via a Hopf bifurcation, the stable manifold disappears, the heteroclinic connection is destroyed, and the corresponding analytical travelling-wave solutions lose their dynamical validity within the full phase-space dynamics. In physical terms, this corresponds to the breakdown of steady wildfire propagation and the transition toward oscillatory, unstable, or spontaneously igniting combustion regimes.

5. Conclusions

In this study, we extended a hyperbolic reaction-diffusion model of a forest fire with relaxation effects by including an advective transport term that accounts for the wind-induced propagation of the fire wave. In fact, the inclusion of an advective term in the model is not just a technical extension, but a qualitative change in the framework of the modeled physical process. In the original reaction-diffusion model, the transport is isotropic and the dynamics of the traveling wave front is controlled primarily by the competition between reaction and diffusion at a finite velocity (via the relaxation term). Once advection is included, the transport becomes directional and the dynamics acquire an explicit dependence on the relative orientation between the wind field and the spreading front. In our extended model formulation, this directional dependence is encapsulated by the oriented effective velocity $v_{\mathrm{eff}}=\mathbf{v}·\ell$, which acts as a key bifurcation parameter controlling both the existence and the qualitative structure of traveling wave fire fronts. Thus, advection alters the propagation regimes of the model by introducing an inherent “preferred direction” and by shifting the stability limits and allowable wave profiles in the parameter space. The main assumption regarding the extended wildfire model presented in this work is that the wind can be oriented either in the direction of the front propagation (downwind) or against it (upwind). This is important because real fire propagation can occur in complex wind-front configurations, including opposing flows or changing wind direction relative to the front normal.

Using SEsM (1,1) with a simple Riccati-type equation, we derived several exact solutions of the model equation that describe the transitions between fire-free and fully fire-affected states or vice versa. Expressing the coefficients in the exact solutions of the model in terms of the coefficients of the simple equation used allowed a unified representation of several physically relevant front profiles within a single analytical framework. Through numerical comparison, we showed that exact solutions with different analytical forms can provide additional information about the thickness, steepness, and spatial attenuation of the fire front and thus serve as comparative configurations that can help further calibrate the model and analyze its sensitivity. Numerical simulations of one of the obtained exact solutions confirmed that variations in the effective wind velocity significantly affect the internal morphology of the fire front.

Furthermore, the analysis of the equilibrium states of the model further elucidates the dynamical mechanisms underlying this process. In a traveling wave system, the existence of an invasion front requires a heteroclinic relationship between the equilibria, which in turn requires one equilibrium to be a saddle while the other is asymptotically stable. Stable wildfire propagation corresponds to the configuration in which the fire-free state $E_0$ is a saddle and the fully fire-affected state $E_1$ is an attractor. Conversely, when $E_0$ is stable and $E_1$ loses stability, sustained fire penetration is dynamically suppressed. In addition, Hopf thresholds, identified by the effective wind velocity, mark critical transitions where the local character of the equilibria changes, altering the geometry of the invariant manifolds and, in certain regimes, destroying the heteroclinic structure that supports stable traveling wave fronts.

The regime map in the $(c,v_{\mathrm{eff}})$–plane can be interpreted as a global framework for different scenarios. On the one hand, it summarizes in a unified geometric way the full range of dynamic possibilities associated with the signs (directions) of c and $v_{\mathrm{eff}}$ , including downwind and upwind propagation, and on the other hand it also accounts for transitions between monotonic and oscillatory regimes and parametric regions where traveling wave fire fronts are not supported. In this sense, the regime map provides a structured basis for exploring “what if” scenarios of fire spread at different wind intensities and orientations, without prior commitment to a single qualitative behavior. In this framework, fire controllability can be interpreted in terms of equilibrium stability: maintaining parameters in regions where the fire-free equilibrium remains stable and no stable heteroclinic relationship exists corresponds to suppression of sustained spread, while regimes in which the fire-affected equilibrium is stable indicate self-sustaining and potentially uncontrolled spread.

At the same time, we emphasize again that the exact traveling wave solutions of the model presented in Section 3 impose an additional analytical constraint connecting the wave speed and the effective wind velocity. In particular, for the families of analytical solutions presented here, there is an additional constraint $c{v}_{\mathrm{eff}}>0$, which means that the projection of the wind onto the front normal and the motion of the traveling wave are co-directed. Therefore, although the global regime map is two-sided and supports both upwind and downwind configurations, the detailed numerical illustrations in this study are limited to the case where the wind and the wave front are moving in the same direction. This constraint should not be interpreted as an inherent limitation of the dynamical system itself, but rather as a characteristic of the particular families of analytical solutions that we have chosen to show in this paper. However, another choice of solution to the algebraic system (21) or another choice of simple equations of the type (10) can allow for the reverse orientation $c{v}_{\mathrm{eff}}<0$, which is a natural direction for future research.

Overall, the exact traveling wave solutions of the extended model proposed in this paper, combined with stability analysis of its equilibrium states and bifurcation analysis, aim to emphasize that wind-driven advection introduces qualitative dynamical transitions, rather than simply quantitative adjustments to the propagation speed of fire—affected tree density. The effective wind velocity emerges as a key control parameter governing the morphology, stability, and admissibility of fire-front waves.

Future work may extend the present framework in several directions. First, an analytical treatment of regimes with opposite wind–front orientation ($c{v}_{\mathrm{eff}}<0$) would allow a systematic investigation of their structure and stability, as well as their implications for wildfire dynamics. Second, incorporating spatially and temporally varying wind fields through the addition of a coupled evolution equation would enable the model to account for dynamically changing environmental conditions and thereby yield a more realistic representation of wildfire propagation.