Efficient Wildland Fire Simulation via Nonlinear Model Order Reduction

We propose a new hyper-reduction method for a recently introduced nonlinear model reduction framework based on dynamically transformed basis functions and especially well-suited for advection-dominated systems. Furthermore, we discuss applying this new method to a wildland fire model whose dynamics feature traveling combustion waves and local ignition and is thus challenging for classical model reduction schemes based on linear subspaces. The new hyper-reduction framework allows us to construct parameter-dependent reduced-order models (ROMs) with efficient offline/online decomposition. The numerical experiments demonstrate that the ROMs obtained by the novel method outperform those obtained by a classical approach using the proper orthogonal decomposition and the discrete empirical interpolation method in terms of run time and accuracy.


Introduction
Accurate modeling of the behavior of wildland fires is an important but challenging task. The models typically involve different physical domains, see [32,37] for an overview, or in parts even data-integrated models [27]. Besides the multi-physics character of the problem, different scales of the involved quantities further complicate the simulation of the coupled models. Consequently, high-fidelity models are computationally costly, conspiring against quick forecasting results required in emergency situations. To reduce the computational time while retaining the complexity of the model, model order reduction (MOR) is a promising approach that has been applied successfully to many different applications. For an overview of MOR methods and applications, we refer to [1,2,5,21,39,40]. Standard projection-based MOR methods such as proper orthogonal decomposition (POD), balanced truncation, or rational interpolation, perform MOR by choosing suitable low-dimensional subspaces of the solution space and then project the system dynamics onto these subspaces using a Petrov-Galerkin approach. The best subspace of a given dimension is classified by the Kolmogorov n-widths [24,38] (respectively the Hankel singular values [50]), thus providing a lower bound on the approximation error. Although in many applications the Kolmogorov n-widths decay exponentially [30], thus enabling MOR to succeed, this is not true for every application. A prominent example is the linear wave equation [17], whose (dimensionless) snapshot matrix based on a Gaussian initial condition is presented in Figure 1a. Comparing this snapshot matrix with the solution of a wildland fire model that   Figure 1b) with the available supply mass fraction (cf. Figure 1c) suggests that the n-widths do not decay exponentially for this model either.
In particular, we cannot expect standard MOR methods to produce an accurate and, at the same time, low-dimensional approximation for wildland fire simulations. To overcome slowly decaying Kolmogorov n-widths, MOR methods for transport-dominated phenomena have emerged over the last couple of years. Let us refer to [7, 9, 25, 28, 34-36, 44, 48, 49, 51] and the references therein to name just a few. As a general strategy, many of these methods use some kind of nonlinear projection framework to allow the ansatz space to adapt along with the solution.
In this paper, we apply the method presented in [7], which can be understood as a generalization of the moving finite element method [33] to general ansatz functions, to a simple wildland fire model introduced in [32]. This showcases the potential of MOR methods tailored to transport phenomena in real-time wildland fire simulations. Our framework builds on an efficient (offline) decomposition [41,42,47] of already computed or measured solution data and constructs the reduced-order model (ROM) via residual minimization as detailed in [7,28]. Due to the nonlinearity of the projection method, the resulting ROM still depends on the original dimension, which conspires against an efficient simulation. Besides, the wildland fire model itself is nonlinear due to a (modified) Arrhenius law, which provides an additional challenge for MOR, even if standard methods are applied. To avoid scaling with the full dimension, it is common practice to work with a further approximation of the nonlinearity, typically referred to as hyper-reduction. A prominent hyper-reduction methods are the empirical interpolation method (EIM) [4] or its discrete analog (DEIM) [12]. In this paper, we extend this approach to cope with dynamically transformed modes. Our main results are the following: • In section 4.2, we extend DEIM by approximating the nonlinearity of the full-order model (FOM) by a linear combination of dynamically transformed ansatz functions to account for the advective transport within the system. Furthermore, we discuss how the state-dependent ROM coefficient matrices can be replaced with cheap to evaluate approximants, see section 4.1. Altogether, the proposed methodology allows for constructing parameter-dependent ROMs while achieving an efficient offline/online decomposition. • Based on the new approach, we construct low-dimensional parameter-dependent ROMs for a wildland fire model, see section 5. Depending on the initial condition, the model inherits complex dynamics that do not allow for a simple separation of transported and non-transported effects, thus rendering this a challenging benchmark problem. Following ideas from [14], we propose a switching strategy, using our method solely in the transport-dominated regime, see section 5.3 for further details. The constructed ROMs allow for accurate predictions within the parameter space and prove to be faster and more accurate than POD-DEIM ROMs, which are based on classical linear subspace approximations. To render the manuscript self-contained, we first present the wildfire model from [32] in section 2. Classical projection-based MOR for a general nonlinear dynamical system is reviewed in section 3.1. Hyper-reduction via DEIM is briefly summarized in section 3.2. In section 3.3, we recall the shifted POD (sPOD) [42] as a powerful offline decomposition method together with the construction of a corresponding ROM based on dynamically transformed modes [7].
Notation. The space of real m × n matrices is denoted by R m×n , and for the transpose of a matrix A, we write A . Furthermore, we use the symbol I n for the identity matrix of size n and denote its entries by the Kronecker delta δ ij = [I n ] ij for i, j = 1, . . . , n. For diagonal matrices and blockdiagonal matrices, we use the abbreviations diag(a 1 , . . . , a n ) : respectively, where a 1 , . . . , a n denote scalars and A 1 , . . . , A n matrices of arbitrary size. In addition, the Kronecker product between two matrices A and B is denoted by A ⊗ B. The symbols ·, · and · are used for the Euclidean inner product and norm, respectively. Furthermore, the space of square-integrable functions over the interval (a, b) ⊂ R is denoted by L 2 (a, b).

Problem Setting -Wildfire Model
We use the model from [32]. Considering a Lipschitz domain Ω ∈ R d with d ∈ {1, 2}, a finite time horizon 0 < t f < ∞, and setting T := [0, t f ], one wants to determine the temperature T : T × Ω → R and fuel supply mass fraction S : T × Ω → R satisfying the coupled partial differential equation (PDE) with reaction rate constant γ S r(T, β, T a ) given by a modified Arrhenius law of the form Here, k denotes the thermal diffusivity, α a proportionality constant representing the ratio of temperature rise per second and reaction rate, β the proportionality coefficient in the modified Arrhenius law, γ a scaled heat transfer coefficient for the heat loss to the environment, γ S the pre-exponential factor, T a the ambient temperature, and v the wind speed (given by atmospheric data or model). An overview of the variables, including the units, is provided in Table 1. For notational convenience, we perform a variable substitution and work with the relative temperature T − T a instead of the temperature, which formally results in T a . For the derivation of the model, fire in a ground layer with finite, small thickness is considered. The fire layer consists of fuel and air above the fuel. The (initially twodimensional) model is then derived from physical principles, namely from conservation of energy in the fire layer and a balance equation for the fuel supply. For a detailed discussion of the derivation, we refer to [32,Section 4]. One of the modeling assumptions is that the only heat source stems from the heat generated by the chemical reaction of burning. Furthermore, it is assumed that heat loss to the atmosphere is only caused by convection and radiation. In contrast, the short-range heat transfer due to radiation and turbulence is modeled via diffusion. More precisely, the individual terms in the model correspond to the following physical effects: • The term ∇ · (k∇T ) is a diffusion term modeling the short-range heat transfer by radiation and turbulence. • The term v · ∇T accounts for the advective heat transfer caused by the wind with wind speed v. • The loss of fuel due to burning is modeled by the product of the supply mass fraction and the reaction rate constant, which is modeled with the modified Arrhenius law (2.2). • The term αγ (T − T a ) models heat lost to the atmosphere due to convection. It is assumed that heat loss to the atmosphere due to radiation is negligible compared to heat loss due to convection, such that the previous term is sufficient to account for both effects.
We remark that the reason for the usage of the modified Arrhenius law lies in the socalled cold boundary difficulty [6]: If a combustion model does not enforce the absence of combustion at ambient temperature, the fuel supply at any given point in space is slowly depleted due to this permanent combustion effect, even before the combustion wave started by some local ignition process reaches this point. Note that the reaction rate constant with modified Arrhenius law is smooth in T . Physical effects that are not accounted for in this model are the effect of fuel disappearance on the total heat capacity of the fuel layer, storage of heat in the ground, chemical kinetics of intermediate products of combustion, the fine scale dynamics of the fire on the surface of and within the wood (fuel), and evaporation of moisture, cf. [31]. The model also does not include any two-way coupling with an atmospheric model. Lastly, note that the Arrhenius law is in general only valid for premixed flames [16] with a sufficient oxygen supply and ignores fuel surface effects.

Preliminaries
3.1. Projection-based MOR and proper orthogonal decomposition. To simplify the following discussion, we assume that we have a finite-dimensional approximation of the coupled PDE (2.1) available, obtained via semi-discretization in space, given bẏ with matrix A T and nonlinear function f . We assume that the approximation error is sufficiently small compared to the model reduction error to come, and refer to (3.1) as the truth model. Combining the discretized temperature T h and supply mass fraction S h in a single state variable z, and the remaining variables in the parameter vector µ, we are thus faced with the task of determining z : T × M → R n for some µ ∈ M with compact parameter set M satisfying Depending on the spatial discretization method, the models in (3.1) and (3.2) may have symmetric positive definite mass matrices in front of the time-derivatives. To simplify notation, we decided to ignore the mass matrices and emphasize that the discussion to follow can easily be extended to the case where the mass matrix is not equal to the identity matrix, by using a suitable weighted inner product. As a further simplification and in agreement with the wildland fire model (2.1), we assume that the matrix A(µ) has an affine decomposition with respect to the parameter, i.e., for some k ∈ N and ν = 1, . . . , k there exist real functions q ν : M → R and matrices A ν ∈ R n×n satisfying The task of projection-based MOR is to identify a suitable low-dimensional subspace Z ⊆ R n , parametrized via an orthonormal basis (ϕ 1 , . . . , ϕ r ) ∈ Z r , such that the solution of (3.2) approximately evolves within Z for any µ ∈ M. Assembling the matrix V = ϕ 1 . . . ϕ r ∈ R n×r and performing a Galerkin projection, we define with projected initial value z 0 := V T z 0 . The approximation obtained by the ROM is then given as z(t; µ) ≈ z(t; µ) := V z(t; µ). (3.5) The question that remains to be answered is how to determine the approximation space Z, respectively the matrix V . If solution trajectories z(·; µ σ ) are available for σ = 1, . . . , with parameter values µ σ ∈ M, then the proper orthogonal decomposition (POD), see for instance [19], constructs the reduced basis by solving the minimization problem (3.6) The POD basis can be computed efficiently via a truncated singular value decomposition (SVD) of the snapshot matrix, given by concatenating the solution vectors at different time instances, which are typically given from a time discretization. We emphasize that the approximation quality throughout the parameter set M will depend on the choice of the training parameters µ σ (σ = 1, . . . , ). If an error estimator is available, then the parameters may be determined iteratively by the POD-greedy method, see [20] and the references therein.
Remark 3.1. MOR for coupled systems of partial or ordinary differential equations is often treated with separate basis functions for the different physical variables, as for instance, the temperature T and the supply mass fraction S in (2.1). This idea of splitting the basis has been utilized in various applications, see for instance [10,22,23,29,52]. Compared to an approach that does not explicitly distinguish between different physical variables, this splitting is often observed to be advantageous in terms of preserving system properties such as stability and passivity.

Efficient approximation of nonlinear terms.
To ensure the usability of the ROM (3.4), one is interested in an efficient offline/online decomposition. This means that in the first stage, the offline phase, one can spend some time constructing the ROM (3.4). In contrast, in the second stage, the online phase, only a little time is available to compute an (approximate) solution for many different parameter values via the ROM (3.4). Analyzing the ROM (3.4), we observe that the matrices A ν ∈ R r×r (ν = 1, . . . , k) can be computed as soon as the reduced basis is determined. Hence, the matrices can be constructed in the offline phase. The evaluation of the linear term A(µ) z thus scales with the dimension of the reduced basis and can be computed efficiently in the online phase. The situation is different for the nonlinear part, encoded in F . Although formally, this is a mapping from R r × M to R r , its computation requires the evaluation of F, thus scaling with the dimension n of the original model, the so-called lifting bottleneck. To avoid this issue, different approaches are proposed in the literature, among them [3,4,11,12,18,26,43].
Here, we focus on the DEIM [12,15], which is the discrete counterpart of EIM [4]. For a given matrix U ∈ R n×m of rank m, the DEIM approximation of F (z, µ) (cf. [12,Def. 3.1]) on span(U ) is given by where S ∈ R n×m is a matrix obtained by selecting certain columns of the n × n identity matrix. The matrix U is typically chosen by finding a best approximation of F (z, µ σ ) for σ = 1, . . . , , i.e., by solving the POD optimization problem (respectively computing the SVD) for the data F (z, µ σ ). Note that the DEIM approximation indeed corresponds to an interpolation along the rows selected by S , since S F DEIM (z, µ) = S F (z, µ). Within the ROM, the approximation of the nonlinear term is thus given by Let us emphasize that the computation of F DEIM requires only certain rows of F , such that an efficient numerical implementation of F DEIM that is independent of the full dimension is possible, thus avoiding the lifting bottleneck. For details we refer to [12]. The question that remains to be answered is the choice of the selection matrix S. While traditionally the columns of I n are selected via a greedy search [4,12], it was recently discovered that it is favorable to compute it by applying a QR-decomposition with column pivoting to U , see [15] for further details.

Nonlinear MOR via transformation operators.
The main drawback of the MOR approach based on linear subspaces is its intrinsic separation of space and time, which cannot always be justified in applications. Prominent examples with strong space-time coupling are the advection equation and the wave equation. Given the snapshot matrices presented in Figure 1, separation of space and time implies that structures along the two principal directions can be approximated efficiently, while diagonally evolving structures are difficult to approximate by a small number of dyadic products. To remedy this issue, one may want to allow the modes ϕ i to evolve over time along with the solution. In [7], this is achieved by replacing the approximation (3.5) with with transformation matrices T i (p i ) ∈ R n×n and unknown path variables p i . For the sake of presentation, we assume p i (t; µ) ∈ R. The optimization problem (3.6) with the modified ansatz (3.9) is then given by i.e., multiple modes ϕ ρ,i are transformed by the same matrix T ρ (p ρ ). Although first numerical studies suggest that it is also possible to optimize over the transformation matrices T i , see for instance [25,45], we focus here on the special case that these matrix functions are given a-priori. For the wildland fire application discussed above, we will use discrete approximations of the shift operator.

Example 3.2.
To demonstrate the benefits of using transformed modes, we consider the linear wave equation with periodic boundary conditions, whose solution is presented in Figure 1a. Using POD, at least r = 49 modes are required to represent the solution with a relative error of less than 1 × 10 −5 . In contrast, using the transformed modes ansatz in (3.9) with a shift operator as transformation, the solution can be exactly, i.e., with no approximation error, represented with r = 2 (respectively ρ = 2, r 1 = 1, r 2 = 1 in (3.11)), see [7,Ex. 4.4].
If we substitute the ansatz (3.9) we immediately observe that T i (p i )ϕ i has to be differentiable with respect to p i , which we assume in the following. Let us denote this derivative with T i (p i )ϕ i and define the path-dependent matrices for the approximation ansatz (3.9) and for (3.11). Similarly as in the Galerkin framework, a ROM can be obtained by projecting the FOM (3.2) with V (p). To obtain equations for the path variables, we perform an additional Petrov-Galerkin projection with V (p) and when using the approximation ansatz (3.9), or defined as D( z) := blkdiag( z 1 , . . . , z q ) ∈ R r×q with z ρ := z ρ,1 · · · z ρ,rρ for ρ = 1, . . . , q when using the approximation ansatz (3.11). These projection matrices naturally arise when enforcing minimization of the residual, cf. [7, sec. 5]. The corresponding ROM is given as Note that (3.12) is a nonlinear system of equations even if the original model (3.2) is linear. Let us emphasize, that in addition to a reduction in space, the ROM (3.12) typically needs considerably less time steps in the temporal discretization for an accurate solution, see [7, sec. 7.3]. [47]. This can be achieved by adding an additional transformation matrix that equals the identity matrix and setting the corresponding path variable to zero. Note that in this case, no further equation is required for the POD path variable.

Efficient Offline/Online Decomposition
A quick analysis of the ROM based on transformed modes (3.12) reveals that an efficient offline/online decomposition as described in section 3.2 for standard projection-based MOR is not readily available, even if the original dynamics are linear. This is because the transformation operators applied to the modes have to be evaluated for different values of the path variables. Thus, in this section, we discuss further approximation steps that have to be implemented to ensure an efficient offline/online decomposition.
for some matrix P ∈ R n×n , and if the right-hand side of (3.2) is equivariant with respect to this transformation, then the explicit dependency of the ROM (3.12) on the transformation matrix and thus on the path variable can be avoided, see [7,Rem. 6.3] for a similar observation. In this case, an efficient offline/online decomposition can be obtained using the methods discussed in section 3.2. An example for a unitary transformation that satisfies an infinite-dimensional analog of (4.1) is the shift operator with periodic boundary conditions, where the spatial differentiation operator assumes the role of P , cf. [7,Ex. 5.12]. If q > 1, then one may still be able to show that the evaluation of (3.12) can be computed without repeated evaluations of the transformation operators, see [8] for details. This is, however, restricted to exceptional cases and not possible in general.
4.1. Efficient approximation of path-dependent matrices. One significant difficulty in evaluating the ROM (3.12) compared to classical linear model reduction methods is that the matrices V and W , and thus also the reduced matrices, depend on the path variable p, which is one of the ROM state variables. This is problematic since assembling the reduced matrices requires computations scaling with the FOM dimension, which needs to be avoided in the online phase.
To this end, we propose to precompute the reduced matrices depending on the path variables for different values of p and construct approximants for these matrix-valued functions based on the sampled data in the offline phase. In the online phase, we can then replace M 1 , N , M 2 , A 1 , and A 2 by their respective approximants to avoid computations scaling with the FOM dimension. The approximants can, for instance, be constructed via interpolation or via regression. Recall that we assume the affine parameter decomposition (3.3) for the matrix A, leading to the affine decompositions for the reduced matrices. In particular, the strategy outlined above is applied to the matrices A i,ν for i = 1, 2 and ν = 1, . . . , k.
Concerning the sampling procedure, we make the following observation. Each entry of the reduced matrices depends at most on two components of the path variable p(t; µ) ∈ R q . If q ≤ 2, then we can simply sample the entire matrix in R q . If, on the other hand, q > 2, it might be advantageous to instead sample each entry of each of the matrices separately, thus alleviating the curse of dimensionality for larger q. A further reduction of the sampling space may be achieved via the method of active subspaces [13], see section 5.1 for further details.
For some choices of the transformation T , the range of p values is naturally bounded. This allows us to restrict the sampling to some finite domain. For instance, when considering the shift operator on a finite computational domain Ω ⊆ R with periodic boundary conditions, each entry of p can be restricted to a domain of the same length as Ω.

Remark 4.2.
A notable special case arises when a single transformation T := T 1 = . . . = T q is used, which additionally satisfies (4.1) as well as for some mapping τ : R × R → R. In this case, it is sufficient to sample the matrices on the left-hand side of (3.12) in a one-dimensional space. This is, for instance, the case for the shift operator on a finite domain Ω ⊆ R with periodic boundary conditions. In that case, the transformation operator satisfies

DEIM with transformation operators.
The methodology outlined in the last subsection allows us to achieve an efficient offline/online decomposition for those ROM terms originating from linear terms in the FOM (3.2). So far, we only addressed those nonlinearities that resulted from the transformation operators that render the considered approximation ansatz (3.9) or (3.11) nonlinear. In this subsection, we focus on the second source of nonlinearity: the FOM nonlinearity F , which occurs in the definitions of F 1 and F 2 and prevents an efficient offline/online decomposition of those terms. To this end, we aim to extend the idea of EIM/DEIM to our setting. For notational convenience, we restrict ourselves to the approximation ansatz (3.11), noting that (3.9) can be handled analogously. The first step of DEIM, see the discussion in section 3.2, consists of finding a suitable linear subspace for approximating the FOM nonlinearity F . As outlined earlier, one popular way is to apply POD to snapshots of F . However, when considering systems with transport, this transport does not only affect the snapshots of the state z, but often also the snapshots of the nonlinearity (see for instance Figure 4 and [46]). As a consequence, we propose a similar approximation ansatz for the nonlinearity as for the FOM state, i.e., we search for ψ ρ,i ∈ R n and b ρ,i : T × M → R for i = 1, . . . , m ρ and ρ = 1, . . . , q with m := q ρ=1 m ρ n and with U (p) := T 1 (p 1 )ψ 1,1 · · · T 1 (p 1 )ψ 1,m 1 · · · T q (p q )ψ q,mq ∈ R n×m , Here, we assume that the transformation matrices and paths used for approximating the state, i.e., T ρ (p ρ ), are also reasonable transformations for approximating the nonlinearity. Especially for the wildland fire model, we expect this assumption to be met, since the nonlinearity is local in the sense that the nonlinearity at a certain point in the spatial domain Ω only depends on the state at the same point. In general, there may be nonlinearities where this assumption is not satisfied and where other paths or even other transformation operators have to be used for approximating the nonlinearity. However, such considerations are not within the scope of this paper. The ansatz functions ψ can be determined based on snapshots of F in the same manner as we determine the ansatz functions ϕ based on snapshots of the state z by (approximately) solving the minimization problem (3.10). Following the EIM/DEIM methodology, we enforce equality in (4.3a) in m ∈ N selected rows. The selected rows within the DEIM framework are typically chosen based on the leading matrix on the right-hand side of (4.3a), i.e., on the matrix U (p). For different values of p, we may thus obtain other row selections, yielding a path-dependent selector matrix S(p). Let us point out that in the online phase, we never explicitly need to assemble the selector matrix itself, but we only need the indices of its non-zero rows. The DEIM approximation with transformed modes, referred to as shifted DEIM (sDEIM), is thus given by As before, we note that, as a consequence of the factor S(p) on the right-hand side of (4.4), we only have to evaluate m rows of the nonlinearity F . Thus, if the Jacobian of the nonlinearity F is sparse, this means that we do not need the full vector V (p) z ∈ R n , but only a small numberm(p) n of its entries. Exploiting this sparsity, we can rewrite the right-hand side of (4.4) as S(p) F (z, µ) = F p, S(p) z, µ . (4.5) Here, F (p, ·, µ) : Rm (p) → R m is the restriction of F to those rows of F required to evaluate the right-hand side of (4.4) and taking only a reduced number of input arguments based on the sparsity pattern of the Jacobian of F . In particular, the function F (p, ·, µ) can be implemented such that its evaluation does not scale with the full dimension n. The details depend on the specific nonlinearity F . They are discussed in section 5.1 for the nonlinearity of the wildland fire model. Within the ROM (3.12), we thus obtain the approximations Since these matrices depend on the path p, we propose, in the spirit of section 4.1, to first compute these matrices in the offline phase for different values of the path p, using, for instance, the QDEIM algorithm from [15]. Afterward, we construct a corresponding interpolant depending on p, which we can then efficiently evaluate in the online phase. At this point, we emphasize that in contrast to the construction of interpolants as outlined in section 4.1, where we never have to sample a space whose dimension is higher than two, here we have to sample a space whose dimension equals q, which is the number of entries in the path p. Again, the dimension may be reduced by exploiting active subspaces [13]. . We refer to [45,46] for similar ideas. Although in general it is possible to extend this idea to multiple transformations, i.e., q ≥ 2, we emphasize that this may result in overlapping DEIM points, which in turn renders the matrix S(p) U (p) singular. This situation did indeed occur in our numerical experiments for the wildland fire, which is the reason we do not pursue this strategy in the following.

Numerical Experiments
In this section, we present the application of the nonlinear model reduction scheme introduced in sections 3.3 and 4 to the wildland fire model (3.1), simulated on a onedimensional strip of length = 1000m, yielding the computational domain Ω = [0, 1000]. To this end, we first discuss how the general nonlinear model reduction framework is applied explicitly to the wildland fire model in section 5.1. Afterward, we present numerical results for two different initial values in sections 5.2 and 5.3. The initial value in section 5.2 leads to a pure advective transport of the already formed combustion waves, whereas the case considered in section 5.3 also includes their formation. In more detail, we use  Table 2 (taken from [32, sec. 7.1]), and initial condition (5.1) for 700s. The solution at the final time is then taken as initial value for our numerical study in section 5.2. The time horizon is chosen so that in all numerical simulations to follow, the combustion waves never reach the boundary of the computational domain. For simplicity, we thus use periodic boundary conditions. Furthermore, as discussed in section 2, we use the relative temperature T − T a for the implementation instead of the absolute temperature T . For the spatial discretization, we use the finite difference method with a 6th order central finite difference stencil. To this end, we decompose the spatial domain into n x = 3000 equidistant intervals corresponding to a mesh width of 1 3 m, yielding a total of n := 2n x = 6000 degrees of freedom (DOF). The spatially discretized nonlinearity in (3.1) is given by where denotes the Hadamard product andr is defined as an entry-wise application of r, see (2.2) The time integration is carried out using MATLAB's ode45 function, which is an explicit time integration scheme with adaptive time stepping based on the Dormand Prince method. For the final time, we use t f = 1400s in section 5.2 and t f = 2100s in section 5.3.
In the numerical experiments, we construct parametric ROMs, which allow variation of the Arrhenius coefficient β. To this end, we use snapshots based on the values β ∈ {540K, 560K, 580K} in the offline phase and test the ROM for various β values from (540-580) K. For notational convenience, we refer to the POD ROM (3.4), where the nonlinearity is approximated as outlined in section 3.2 via DEIM, as POD-DEIM ROM. Accordingly, our method resulting in the ROM (3.12) with approximation of the nonlinear terms as in section 4 is abbreviated with sPOD-sDEIM ROM.
All relative errors stated in the following subsections are measured in a discretized L 2 (0, t f ; R n ) norm, where we used the trapezoidal rule for approximating the integral. 5.1. MOR for wildfire model. In this subsection, we discuss how we apply the general model reduction procedure outlined in sections 3.3 and 4 to the wildland fire model (3.1). Since the (relative) temperature and the supply mass fraction have different scales, we follow the principle described in Remark 3.1 and determine different basis functions for the temperature and for the supply mass fraction, respectively. For each of these two quantities, we use (3.11) as approximation ansatz with q = 2 for the two traveling combustion waves as shown in Figure 1b and Figure 1c. Since the propagation of the combustion waves is reflected in both the snapshots of the temperature and the snapshots of the supply mass fraction, we take the same paths for the temperature and supply mass fraction, i.e., in total, we need to compute two paths in the online phase. The method we use to compute the basis functions according to the approximation ansatz (3.11) is outlined in sections 5.2.1 and 5.3.1, respectively. For transforming the basis functions, we use the constant extrapolation shift operator, which has been applied in [47] for constructing lowdimensional descriptions of reaction fronts in the context of combustion. Let us emphasize that the supply mass fraction attains different (approximately) constant values before and after the reaction front, such that a shift operator with periodic boundary conditions would lead to unphysical supply mass fraction modes, which do not properly capture the shape of the respective reaction front.
Whenever the shift operator needs to be evaluated at points that are not multiples of the spatial mesh width, we use polynomial interpolation based on a cubic Lagrange polynomial as in [7, sec. 7]. This discretization of the shift operator leads in general to a non-differentiability of the mapping p → T (p)ϕ at multiples of the spatial grid size. To remedy this issue, we propose to first compute the derivative of this mapping on the infinitedimensional level and discretize afterward. As a result, we obtain T (p)ϕ := T 0 (p)Dϕ, where T 0 (p) is a discretization of the shift operator with zero extrapolation and D is a sixth-order central finite difference approximation of the first derivative with one-sided finite differences at the boundary.
In the offline phase, the paths are determined by tracking the highest positive and negative value of the discrete spatial derivative within the temperature snapshots. Here, the highest positive value corresponds to the front of the left-going combustion wave and the highest negative value to the right-going one, cf. Figure 12 for a depiction of the combustion waves. Since this approach yields a step function for the paths, we partially replace the path via a linear interpolation in a post-processing step to obtain smooth functions. For the initial condition considered in section 5.2, the complete path trajectories are replaced by linear interpolants, which are based on the first and on the last value of the path with respect to time. This reflects that the propagation speeds of the combustion waves are approximately constant. The case considered in section 5.3 additionally covers the formation of the wave profiles. Therefore we replace the path within the last 98.5% of the time interval by a linear interpolant. In contrast, the empirically determined path trajectory in the first 1.5% of the time interval is not a linear function of time. It is thus taken as obtained via the front-tracking approach.
To be able to evaluate the terms M 1 , M 2 , N , A 1 , and A 2 occurring in the ROM (3.12), as well as V , W , and V in (4.6) efficiently, we use the approach presented in section 4.1 as follows: In a first step, we compute an SVD of the path variables determined in the offline phase to determine an active subspace. We observe that it is sufficient to sample within the subspace showing that the two combustion waves propagate with the same velocity in opposite directions. This is due to the symmetric initial condition. Exploiting this symmetry, we only need to sample a one-dimensional domain [0, ], and based on the sampled data, we construct a piecewise constant interpolant in the offline phase, which can be evaluated efficiently in the online phase. The nonlinearity of (3.1) is given by where ⊗ denotes the Kronecker product. Since the nonlinearity is given by a Kronecker product of a constant vector and a nonlinear function, we only need to approximate f (S h , T h , β) instead of F (z, µ) by the sDEIM method outlined in section 4.2. To this end, we use again the approximation ansatz (3.11) for determining a low-dimensional description for f (S h , T h , β), which yields U (p) in the notation of section 4.2. The rest follows along the lines of section 4.2, with the only difference that the products V (p) U (p) and W (p) U (p) in (4.6) have to be replaced by respectively. The incidence matrix corresponding to the sparsity pattern of the nonlinearity f (S h , T h , β) is given by i.e., the ith entry of f depends only on the ith and on the (i + n x )th entry of z for i ∈ {1, . . . , n x }. As a consequence, the number of required rows of V (p) in (4.5) ism = 2m and, thus, not dependent on the path. Hence, we can achieve an efficient offline/online decomposition of the term S(p) V (p) as outlined at the end of section 4.2. Finally, the modified nonlinearity F from (4.5) is here simply given by T 1 · · ·T mŜ1 · · ·Ŝ m → Ŝ 1r (T 1 , β) · · ·Ŝ mr (T m , β) .  Especially, we note that for the specific nonlinearity of the wildland fire model, F does not depend on the path. Moreover, the mapping F allows to vary the parameter β in the ROM without requiring a parameter-affine structure in f .

Initial condition with traveling wave solution.
In this section, we discuss the numerical solution of the ROM (3.12) for the initial condition where the initial ignition process has already taken place, and combustion waves have already developed.

Determination of modes.
We use the following procedure to approximately solve the optimization problem (3.10): We first separate the two combustion waves by dividing the snapshot matrix in the middle (cf. Figure 2a) and replacing the left/right combustion wave with zeros, see Figure 2b. Afterward, each set of snapshot data is now transformed to the co-moving frame, which is obtained by shifting the spatial domain according to the pre-determined path. The resulting transformed snapshot data is low-rank, and suitable modes are now determined from the transformed snapshot data via POD, respectively SVD. The approximations for each co-moving frame are then shifted back and added back together, to get an approximation of the combined snapshot data. The corresponding offline errors and online errors for Arrhenius coefficient β = 558.49K for different numbers of modes are presented in Table 3, showing that this strategy seems indeed reasonable to tackle the optimization problem (3.10). Note that the DOF are given by the total number of transformed modes plus the additional two path variables, which we indicate with the notation in Table 3. The first two modes for the left-going combustion wave are presented in Figure 3. We see, in agreement with the offline errors presented in Table 3, that already the first mode captures the dominant dynamics. Let us emphasize, that the coefficient for the first temperature mode varies less than 2 × 10 −3 percent, while the coefficient for the first supply mass fraction mode varies less than 9 × 10 −3 percent.

Efficient offline/online decomposition.
For the sDEIM method as described in section 4.2, we first notice that the snapshots of the nonlinearity indeed feature similar combustion waves with the same wave speeds, see Figure 4 for an illustration. This justifies that we take the same transformations and the same paths used for the approximation of the state to approximate the nonlinearity. Moreover, the mode numbers are chosen  such that the number of nonlinearity modes is twice the number of modes used for the temperature and supply mass fraction (both temperature and supply mass fraction have the same number of modes). As outlined in sections 4.1 and 4.2, we need to sample a couple of path-dependent matrix functions to achieve an efficient offline/online decomposition. Therefore, as mentioned before, we make use of the active subspace defined in (5.2), and only sample within a finite subset of this subspace. For the settings considered in this subsection, it appears to be sufficient to use {(p 1 , p 2 ) ∈ U a | p 1 ∈ [0, 300]} as sampling domain, since the combustion waves never travel more than 300 meters for the considered range of Arrhenius coefficients β and for the considered time horizon. For the discretization of the sampling domain, we use an equidistant sampling grid with a grid size of 20 3 meters, since this leads to reasonable approximation errors and computation times for all considered mode numbers.

ROM simulations.
Using the efficient offline/online decomposition, we can now compare the efficiency of the ROMs for different numbers of modes. To this end, we always use the same number of modes per path or reference frame, and we use the same number of modes for the temperature and for the supply mass fraction. The total number of DOF for the ROM obtained by our approach is then given as r + 2, where r is the total number of transformed modes and 2 is the number of frames for the considered test case, i.e., the number of path variables. The details of the errors and speedups for different   Table 3, showing that the relative online errors are within one order of magnitude to the offline errors. We notice that while the error in the offline phase decays with an increasing number of modes, this is not necessarily true in the online phase. Let us emphasize that the online error is subject to several different approximations, such that we cannot provide a decisive explanation for the observed behavior. A possible explanation for stagnation in the online error seems to be that with an increasing number of modes, the ROM is more sensitive to errors stemming from the spatial discretization. A detailed error analysis is not within the scope of this paper and subject to future research.
We compare POD-DEIM ROMs and sPOD-sDEIM ROMs for the Arrhenius coefficient β = 558.49K with respect to computation time and accuracy, for different numbers of modes. To avoid irregularities in the computation time, we average the computation time over three simulations per distinct mode number. The results are presented in Figure 5. We find that the sPOD-sDEIM ROM outperforms the POD-DEIM ROM, both with respect to accuracy and computation time. Note that the POD-DEIM ROM with 200 modes still results in an error that is one order of magnitude worse than all sPOD-sDEIM ROMs.
To investigate the ability of the sPOD-sDEIM ROM to make predictions when evaluated at parameter values not in the training set, we use the sPOD-sDEIM ROM with 3 modes per variable and frame, yielding a total of 3 × 2 × 2 + 2 = 14 DOF (including the path variable for each frame). As training parameter values, we use three samples of the FOM with β ∈ {540K, 560K, 580K}. We test the sPOD-sDEIM ROM for various β values within the range from 540K to 580K. The resulting relative online errors for the temperature and for the supply mass fraction are depicted in Figure 6. For computing the errors, we also simulated the FOM for these parameter values. We observe that with the settings at hand, the ROM is able to make accurate predictions over the complete sampling interval, with worst-case error less than 1%. In average, our computations show that the ROM is more than 75 times faster than the FOM with a mean error of less than 0.5%.

Initial condition with combined traveling and non-transported effects.
We now turn to the more challenging scenario with the initial value given in (  -Relative online error of the sPOD-sDEIM ROM with 12 + 2 DOF for β ∈ {540K, 560K, 580K} over the Arrhenius coefficient β ranging from 540K to 580K and 1c. Let us emphasize that in the beginning of the evolution of the temperature variable (cf. Figure 7a), a separation of the traveling waves and the ignition is not straightforward. Instead, we follow the strategy from [14] and partition the time interval to construct two different surrogate models: • In the beginning of the time interval, i.e., in the area (i) in Figure 7a, we approximate the wildland fire model with standard POD-Galerkin, as described in section 3.2. • As soon as a clear separation of the waves and the ignition is possible, i.e., in areas (ii-a) and (ii-b) in Figure 7b, we use our ROM with transformation operators to capture the two wave fronts and add additional POD modes, as described in Remark 3.3, for the remaining dynamics.

Determination of modes.
As outlined earlier, we use a POD approximation in the time interval t ∈ [0, t switch ], which corresponds to the area (i) in Figure 7b, and a combined approximation of transformed modes and POD modes in the interval t ∈ [t switch , t f ] with t switch = 100s, corresponding to areas (ii-a) and (ii-b).  We have to separate the traveling waves and the non-transported effects for the combined approach with transformed modes and POD modes. To minimize the impact of the nontransported effects on the modes for the traveling waves, we further separate the area after the switch and extract the transformed modes solely from area (ii-b), which consists of 65% of the time domain after the switch for the relative temperature and 80% for the supply mass fraction. We, therefore, propose the following strategy: (i) In area (ii-b) in Figure 7b, we proceed with the heuristic discussed in section 5.2.1, i.e., we separate the two combustion waves by dividing the computational domain in the middle, replacing the missing parts with zeros, and then shift accordingly.
For an illustration, we refer to Figure 2. In the shifted frame, i.e., in Figure 2b, we apply standard POD to obtain the transformed modes. (ii) To eliminate the traveling waves also from the area (ii-a), we fit a low-degree polynomial to the coefficients z ρ,i (t; µ) in (3.11) for the traveling wave modes obtained within the shifted frames, and use these polynomials to extrapolate the coefficients. (iii) We subtract the traveling wave approximations from the snapshot data and apply standard POD to capture the non-transported effects.
We conclude the offline phase with several important remarks: • The above strategy (and similarly the method in section 5.2.1) will, in general, not provide a minimizer for the optimization problem (3.10). Nevertheless, these heuristic methods can be computed efficiently and deliver satisfactory results. In contrast, techniques that are based on solving the minimization problem (3.10) or related minimization problems, cf. [41,47], involve expensive iterative solvers for an optimization problem, whose number of parameters scales with the dimension of the spatial and temporal discretization. • Note that we extrapolate the coefficients instead of computing them via projection.
The main reason for this is that if we simply project, then the non-transported effects are, in parts, also approximated by the transformed modes. This can be We observe different error decays for the temperature and the supply mass fraction: The error decay in the temperature suggests that the numbers of POD modes chosen for the first and second time intervals should be adequately balanced. For instance, if only one POD mode per variable is chosen for the second time interval, it does not pay off to increase the number of used POD modes for the first time interval to values higher than 10, since for higher values, the errors stagnate. In contrast, the error decay for the supply mass fraction seems to indicate that the error is more or less independent from the number of POD modes used in the second time interval. These different behaviors in the temperature and supply mass fraction error are likely because at the beginning of the second time interval, there is still a fading temperature peak in the middle of the computational domain. However, this peak is not reflected in the supply mass fraction, which is already almost entirely consumed in the middle of the computational domain at the beginning of the second time interval. Consequently, adding POD modes in the second time interval is especially valuable for approximating the temperature but less significant for the supply mass fraction. Despite this observation, we decide, for simplicity, to always take the same number of modes for the temperature and for the supply mass fraction while noting that there is some unexploited potential for further improvement. temperature and for the supply mass fraction, whereas we use twice as many modes for approximating the nonlinearity f . For the second time interval, we use a ROM with transformed modes to capture the traveling combustion waves, enriched with additional POD modes for describing the fading combustion in the middle of the computational domain. Again, we use the same number of modes for the temperature and for the supply mass fraction, whereas we choose twice as many POD and transformed modes for the nonlinearity. For sampling the path-dependent matrices, we exploit the active subspace as in section 5.2.2 and sample within the range {(p 1 , p 2 ) ∈ U a | p 1 ∈ [0, 500]} . This sampling domain is discretized using an equidistant grid with a grid size of 1 3 meters. In between the two time intervals, we need to switch between the two ROMs. This requires computing a new initial value for the ROM used in the second time interval. To this end, we first compute the approximation of the full-dimensional state based on the reduced state calculated at the last time point of the first time interval. Afterward, we use this approximation to determine the initial paths of the second time interval by using the same tracking procedure as in the offline phase, see section 5.1. Finally, the remaining coefficients of the initial condition are computed by an orthogonal projection of the approximation onto the span of the transformed modes and the POD modes used for the second time interval.

ROM simulations.
We conclude our numerical study by detailing the online performance indicators of the sPOD-sDEIM ROM. We start by comparing a standard POD-DEIM ROM with our method in terms of computation time and accuracy. The results are presented in Figure 10. As before, the computation times are averages over three simulations. As in the previous section, our method outperforms standard POD with a better accuracy to computation time ratio for all considered DOF.
Similarly as before, we perform a parameter study for the sPOD-sDEIM ROM by changing the Arrhenius coefficient β. To this end, we construct a ROM with 28 POD Arrhenius coefficient β [K] relative online error relative temperature supply mass fraction Figure 11 -Relative online error of the sPOD-sDEIM ROM with 28 POD modes before switching and 12 transformed modes, 4 POD modes, and 2 path variables after switching, trained for β ∈ {540K, 560K, 580K}, over the Arrhenius coefficient β ranging from 540K to 580K modes before the switching time t switch = 100s, using 14 modes for the relative temperature and 14 modes for the supply mass fraction. After the switch, we use 3 transformed modes for each variable and each frame and 2 additional POD modes per variable, giving a total of 18 DOF, including the path variables for each frame. The modes are determined from snapshot data obtained from the simulation of the FOM for the training parameters β ∈ {540K, 560K, 580K}. The relative errors are depicted in Figure 11. Clearly, our ROM can produce accurate results throughout the parameter domain with an average speedup factor higher than 25, see Table 4 for further details. The worst-case approximation error throughout the parameter domain is achieved at β = 540K. The corresponding solution of the FOM and the sPOD-sDEIM ROM for the temperature at different time instances  is depicted in Figure 12, detailing that the wavefronts are captured accurately with a negligible offset.

Conclusions
In this paper, we introduce a hyper-reduction framework for the recently introduced nonlinear model reduction scheme presented in [7] and apply the method to a wildland fire model with nonlinear dynamics, featuring traveling combustion waves. As approximation ansatz, we use a linear combination of transformed modes, cf. (3.9) and (3.11). Thus, the reduced state vector consists of the coefficients of the linear combination and the so-called path variables that parameterize the transformation operators. The hyperreduction framework addresses two sources of nonlinearity in the reduced-order model: First, we consider those nonlinearities originating from the nonlinear approximation ansatz. Primarily, we propose to sample the path-dependent coefficient matrices in the offline phase and replace them with interpolation-based approximants in the online stage. Second, we consider the nonlinearity originating from the FOM nonlinearity. We propose an extension of the discrete empirical interpolation method (DEIM) by approximating the nonlinearity with a linear combination of transformed modes. All in all, this new hyper-reduction framework allows obtaining parameter-dependent ROMs based on transformed modes while achieving an efficient offline/online decomposition.
We apply the new method to two wildland fire test cases with different initial conditions, using a shift operator as a transformation to capture the traveling combustion waves. The first considered initial condition consists of already separated combustion waves, such that the dynamic is mainly characterized by the propagation of the waves. In this setting, the proposed model reduction framework based on transformed modes yields very low-dimensional ROMs, which clearly outperform ROMs constructed by a classical POD-DEIM approach in terms of accuracy and speedup. Furthermore, the ROMs are parameter-dependent and prove to yield accurate predictions in the considered range of parameter values. The second considered test case also includes the formation of the combustion waves, thus combining traveling and non-transported effects. Consequently, this is a challenging problem, not only for standard POD-DEIM ROMs, but also for a ROM based only on transported modes. For this setting, we propose to take a low-dimensional POD-DEIM ROM for the initial time period where the combustion waves develop, and a ROM based on transformed modes and POD modes in the second part of the time interval. Again, the ROM obtained by this approach is significantly faster and more accurate than classical POD-DEIM ROMs. Furthermore, the ROM yields accurate predictions within the considered parameter range.
Despite the convincing observed performance of our approach, there are still some shortcomings and difficulties which need to be addressed in the future. First, we use the shift operator for the presented numerical experiments on a one-dimensional domain, but the extension to higher space dimensions is non-trivial. To be able to cope with complex transports in higher dimensions, one promising direction is to use more involved coordinate transforms, see for instance [25,45,48] for some contributions in this direction. Another difficulty that we observed in the numerical experiments is that the quality of the ROM obtained by our approach relies on a sufficiently fine spatial discretization of the FOM. We assume that the reason for this is that we need to discretize the derivative of the shift operator for constructing the ROM and the corresponding discretization error scales with the spatial discretization error of the FOM, see also [7, sec. 7.2]. Furthermore, in contrast to the classical POD-DEIM approach, our method requires choosing more parameters. Their influence on the error is less evident than for POD-DEIM.
Besides the shortcomings mentioned in the last paragraph, an interesting future research direction is to ensure stability within the sPOD-sDEIM ROM based on transformed modes, for instance, by using a port-Hamiltonian formulation of the governing equations. First investigations suggest that the considered wildland fire model allows for such an energybased formulation, which could be exploited for structure-preserving model reduction schemes. Furthermore, applying the proposed method to more involved wildland fire models is another interesting and challenging task for the future. For instance, we could consider wildland fires formulated on two-dimensional domains or take into account the fluid dynamics in the atmosphere.
funding from the DFG under Germany's Excellence Strategy -EXC 2075 -390740016 and is thankful for support by the Stuttgart Center for Simulation Science (SimTech).