XX October 20XX XX 01 03

# Physics-Informed Machine Learning Simulator for Wildfire Propagation

Luca Bottero*, Francesco Calisto*, Giovanni Graziano*, Valerio Pagliarino*,
Martina Scauda†, Sara Tiengo* and Simone Azeglio* *Student  E-mails:
luca.bottero192, francesco.calisto, giovanni.graziano136, valerio.pagliarino,
Department of Physics  E-mails: luca.bottero192, francesco.calisto,
giovanni.graziano136, valerio.pagliarino,  University of Turin  E-mails:
luca.bottero192, francesco.calisto, giovanni.graziano136, valerio.pagliarino,
†Student  E-mails: luca.bottero192, francesco.calisto, giovanni.graziano136,
valerio.pagliarino,  Department of Mathematics  E-mails: luca.bottero192,
francesco.calisto, giovanni.graziano136, valerio.pagliarino,  University of
Turin  E-mails: luca.bottero192, francesco.calisto, giovanni.graziano136,
valerio.pagliarino,

ABSTRACT

The aim of this work is to evaluate the feasibility of re-implementing some
key parts of the widely used Weather Research and Forecasting WRF-SFIRE
simulator by replacing its core differential equations numerical solvers with
state-of-the-art physics-informed machine learning techniques to solve ODEs
and PDEs, in order to transform it into a real-time simulator for wildfire
spread prediction. The main programming language used is Julia, a compiled
language which offers better perfomance than interpreted ones, providing Just
in Time (JIT) compilation with different optimization levels [1]. Moreover,
Julia is particularly well suited for numerical computation and for the
solution of complex physical models, both considering the syntax and the
presence of some specific libraries such as DifferentialEquations.jl and
ModellingToolkit.jl.

##  1 INTRODUCTION

In recent years wildfires have been increasingly growing in intensity and
frequency, becoming a serious threat to the health and socio-economic
stability of various countries all over the world. In particular, Australia
was devastated by the "Black Summer", the bushfire season between 2019 and
2020, and California just suffered from the most severe wildfire season
recorded in its modern history. According to the California Department of
Forestry and Fire Protection [2]: over 4 percent of its land was burned by
more than 8,600 fires [3].

The interconnection between wildfires and climate change is evident and has
been studied by climate scientists over the years [4][5]. On one hand, climate
change favors the spread of wildfires. One of the main factors in play is the
rise in global temperatures: the early beginning of spring leads to the rapid
melting of snowpacks, causing land to dry out earlier and remain dry for
longer. In addition, the hotter the air, the more water it soaks up from
plants and soils; the "vapor pressure deficit" is used to measure the
difference between how much water the air holds and how much it could hold:
the higher this coefficient, the more soil and vegetation will dry out [6].
Furthermore, bark beetles and other insects that survive at high temperatures
are responsible for killing millions of trees, according to the Fourth
National Climate Assessment [7], turning them into kindling for wildfires.
These two phenomena together produce an alarming growth in wildfire frequency
and intensity or, following the terminology used in [4], the number of fires
(NB) and the extension of the burnt area (BA). On the other hand, wildfires
impact climate change in two ways: firstly, they emit massive amounts of
carbon dioxide and other pollutants that can affect regional and even global
climate. Secondly, vast lands laid bare by fires are soaked with rain,
increasing the risk for devastating landslides. On top of this, wildfires
produce serious health risks due to the emission of toxic chemical compounds
[8].

In such a worrisome scenario, having the support of fast and accurate wildfire
spread simulators that can run iteratively in order to evaluate several
different containment strategies is of unprecedented importance.

###  1.1 Real-time is revolutionary

The current stable version of WRF is able to produce highly accurate and
validated predictions, but each simulation takes several hours, which is
acceptable for a single run. Nevertheless, if one wants to use the simulator
for containment strategies, time is a major obstacle. In order to use WRF for
that purpose, one should run the simulator iteratively with different initial
scenarios, therefore predicting the outcome of various containment strategies
and choosing the best one according to the output of the program. The
revolutionary reach of such an approach is that it would limit human error in
situations where taking the right decision quickly is extremely difficult. To
reach this goal it is essential to significantly increase the computational
efficiency of fire spread models, which can be done thanks to Machine Learning
techniques.

Furthermore, in our proposal we will show that such techniques introduce a
broad spectrum of other advantages, for example: the possibility to obtain a
continuous solution for the problem instead of a discretized one on a grid,
the possibility to forecast the evolution even outside the limits of the
simulation domain, and the possibility to reconstruct the evolution back in
time, with applications for forensic analysis. The main programming language
used in this work is Julia for reasons explained below.

###  1.2 Physics-based or ML?

In the context of science, the well-known adage “a picture is worth a thousand
words” might well be “a model is worth a thousand datasets.”, writes
Christopher Rackauckas in his trademark paper [9]. A single sentence couldn’t
describe the true power of Scientific Machine Learning better. Machine
Learning’s utmost strength, i.e. universal flexibility to approximate any
nonlinearity from data, as stated by the Universal Approximation Theorem [10],
is at the same time its main drawback: the higher the complexity of the
problem, the more data the algorithm needs to train. While some fields can
provide big data, such as bioinformatics [11], many scientific disciplines are
still limited in the use of Deep Learning (DL) by a lack of data. Another
reason for the ineffectiveness of Machine Learning (ML) might be the apparent
chaotic nature of a phenomenon. This is the case of wildfires and weather-
related phenomena, which don’t depend solely on the initial and boundary
conditions, but evolve according to complex internal dynamics that cannot be
predicted from mere data. In order to make a neural network learn such
dynamics, one would need data with both high spatial and temporal resolution,
something not available in this field [12]. For this reason, we think that the
growing field of Scientific Machine Learning is a good compromise where
climate and wildfire science could take advantage. The idea is to embed
physical information, generally in the form of partial differential equations,
into the ML model, thus constructing a Physics Informed Neural Network (PINN).
A former physics-informed approach consists in training the model on data and
constraining the solution space to the physically admissible ones (e.g. for
incompressible fluid dynamics, any flow solution that breaks the mass
conservation principle would not be considered by the model [13]). Another
physics-informed architecture - the one considered in our case - consists
instead in training the neural network directly on the physical model output
data, employing the feed-forward networks as an efficient approximator of
trajectories in the solution space. This approach leads to various advantages,
such as avoiding the so-called curse of dimensionality: as the dimensionality
of an equation increases, the computational costs of a numerical equation
solver grow exponentially, while PINNs can be proven to have a polynomial
bound [14]. This results in a significant speed up with respect to numerical
solvers when solving high dimensional problems. We decided to use the second
technique since our goal is to evaluate the feasibility of a high efficiency,
ML-based, wildfire spread simulator that uses the equations of the WRF-SFIRE
model.

##  2 RESEARCH QUESTION

We aim to investigate the possibility to implement some key parts of the
complex WRF-SFIRE physical model using an ML-based approach with all the
advantages disclosed above and discussed in paragraph 4.4.2.

##  3 METHODS

Firstly, we started by investigating the structure of the systems of partial
differential equations (PDEs) that constitute the core of the physical model
of WRF-SFIRE. In this way we detected the modules where replacing the current
structure has the greatest impact and we found out how to re-implement them
using our architecture.

###  3.1 WRF model exploration and profiling

In order to increase the computational efficiency of WRF-SFIRE model, we
needed to identify the most computationally-intensive parts of WRF’s code. We
performed a profiling using the Linux perf command [15]. In this way, we
obtained the CPU overhead of the subroutines measured running a variety of
notable example scenarios provided by the developers of the WRF model. For
comparison, we also ran a performance analysis for atmospheric simulations
alone. The results we obtained are consistent with our hypothesis about what
slows down the simulations: numerical solvers of the model’s main differential
equations are responsible for more than $20\%$ of the total overhead. The
remaining machine time is due to I/O processes and other equations that
require less complex computations or that are updated less frequently. The
main difference between the purely atmospheric simulations and those with a
fire simulation was, as expected, due to calculations of the fire spread rate
and the resolution of the level set equation, which determines the fire
propagation. This analysis was followed by an in-depth examination of the
WRF’s Fortran code, with particular attention to the SFIRE module, in order to
understand how the problem of solving the governing equations is approached.
We used this information to decide which equations to re-implement in our
model. This is our verdict: the level set equation, describing the evolution
of fire area, is not only central from a conceptual point of view, but its
numerical resolution has a significant impact on the CPU load during fire
simulation. (Please see Fig.1 and Fig.1). We were also interested in
evaluating the possibility of applying PINNs to the resolution of the Euler
system that is a system of 7 PDEs governing the atmosphere behavior. That is
because even if the Euler system is not strictly linked to wildfire
simulation, it has the greatest impact on the computational load and is
physically coupled to the level set equation.

Therefore, in this paper we will present a detailed application of the PINN
architecture to the level set equation and a convergence study of PINNs
training applied to the Euler system.

Figure 1: Profiling of computation-expensive WRF modules Module name (CPU, 4 threads) |  Overhead %
---|---
01\. sfire_phys_MOD_fire_ros_cawfe |  8.70
02\. sfire_core_MOD_check_lfn_tign_ij |  6.66
03\. powf_fma_[Standard library for basic math] |  6.58
04\. small_step_em_MOD_advance |  5.99
05\. sfire_util_MOD_print_3d_stats |  5.89
06\. small_step_em_MOD_advance_uv |  5.29
07\. fs_sfire_phys_MOD_fire_ros |  4.64
08\. small_step_em_MOD_advance_mu_t |  4.48
09\. sfire_core_MOD_tend_ls |  2.55
10\. advect_em_MOD_advect_scalar_pd |  2.52
11\. fr_sfire_phys_MOD_fire_rate_of_spread |  2.41
Table 1:

###  3.2 A brief insight into the WRF-SFIRE model

Despite the focus of this project being on wildfire spread simulation, the
current state-of-the-art fire spread models can produce realistic results only
because they are integrated with a more complete atmospheric simulator;
therefore presenting the structure of WRF-SFIRE is impossible without insight
into the atmospheric model that is interrogated continuously by the fire
propagation module.

###  3.3 Description of the physical model

The Weather Research and Forecasting (WRF) model is an atmospheric modeling
system designed for Numerical Weather Prediction (NWP). The NWP procedure was
designed by the meteorologist Vilherlm Bjerknes in 1904, originating from the
assumption that weather forecasting is an initial value problem. Therefore,
once given the observed initial state of the atmosphere, the integration of
governing equations is feasible and is assigned to the Advanced Research WRF
(ARW), the dynamics solver of the model. However, the equations cannot be
solved analytically, so they are described on a mesh and the computed
solutions are an approximation of exact ones. A good and comprehensive
description of the model is accessible in the Users’ guide [16].

#### Grid and coordinates

The atmospheric model operates on a logically quadrilateral 3D grid on the
Earth’s surface (Fig.2). The vertical coordinate is the hybrid terrain-
following mass coordinate $\eta$; defined as a combination between the
terrain-following coordinate and a pure pressure coordinate. Considering the
hydrostatic pressure $p_{d}$ as an independent variable, ARW integrates the
compressible, non hydrostatic Euler equations cast in the flux form using the
following prognostic variables that have conservation properties ([17], pp.
8-10):

  * •

$\bm{v}=(u,v,w)$, where $\bm{v}$ are the covariant velocities in the
horizontal and vertical directions, while $\omega=\partial_{t}\eta$ is the
contravariant vertical velocity;

  * •

$\theta_{m}$ is the moist potential temperature;

  * •

$q_{m}=(q_{v},q_{c},q_{r}...)$ represents the mixing ratios of moisture
variables (water vapor, cloud water, rain water, …).

Since equations are formulated in flux form, the previous conserved prognostic
variables have to be defined in flux form as well, leading to

$\begin{cases}\bm{V}=\mu_{d}\bm{v}=(U,V,W)\\\ \Omega=\mu_{d}\omega\\\
\Theta_{m}=\mu_{d}\theta_{m}\\\ Q_{m}=\mu_{d}q_{m}\\\ \end{cases}$ (1)

where $\mu_{d}=\partial_{\eta}p_{d}$ is the vertical coordinate metric. On the
other hand, the geopotential $\phi=gz$, despite being a prognostic variable in
the governing equations of the ARW, is not written in flux form because the
quantity $\mu_{d}\phi$ is not conserved.

#### Downscaling and nesting

Even though the aim is to simulate microscale and mesoscale events, the
synoptic and global scales cannot be neglected. However, a high resolution
global model is computationally expensive. Through a downscaling it is
possible to enhance the resolution starting from large scale information. WRF
implementes dynamic downscaling through the nesting method that introduces an
additional grid (child grid) into the simulation with different granularity.
Prognostic variable fields from the coarse-grid forecast are imposed as a
boundary condition of the integration area.

Figure 2: Coordinates system of model’s grid [17]

###  3.4 Fire module

The coupling of a mesoscale weather model (WRF-ARW model) with a 2D semi-
empirical fire spread model permits the execution of a wildland fire
simulation [18]. SFIRE was developed based on the level set method i.e. an
Eulerian approach. The reason why the fire model needs to be coupled with the
atmospheric model is due to the significant influence of weather on fire
behavior. In particular, the wind affects fire spread, while the fire causes
changes in the atmosphere through the heat and vapor fluxes from burning fuel
and evaporation of fuel moisture. For every WRF time step, one time step of
the fire model is performed. The called fire model inputs the surface wind,
which drives the fire, and outputs the heat flux into the atmosphere. The fire
model has a high resolution and operates on a refined fire mesh contained in
the inner grid of the atmospheric model.

###  3.5 Mathematical structure

#### Atmospheric governing equations

Using the variables defined in (1), the flux-form Euler equations,
representing the core of the ARW, can be written as follows [17]

Conservation of momentum

$\begin{cases}\partial_{t}U+(\nabla\cdot\bm{V}u)+\mu_{d}\alpha\partial_{x}p+(\frac{\alpha}{\alpha_{d}})\partial_{\eta}p\partial_{x}\phi=F_{U}\\\
\partial_{t}V+(\nabla\cdot\bm{V}v)+\mu_{d}\alpha\partial_{y}p+(\frac{\alpha}{\alpha_{d}})\partial_{\eta}p\partial_{y}\phi=F_{V}\\\
\partial_{t}W+(\nabla\cdot\bm{V}w)-g[(\frac{\alpha}{\alpha_{d}})\partial_{\eta}p-\mu_{d}]=F_{w}\\\
\end{cases}$ (2)

Conservation of heat

$\partial_{t}\Theta_{m}+(\nabla\cdot\bm{V}\theta_{m})=F_{\Theta_{m}}$ (3)

Conservation of mass

$\partial_{t}\mu_{d}+(\nabla\cdot\bm{V})=0$ (4)

Geopotential Material Derivative

$\partial_{t}\phi+\mu_{d}^{-1}[(\bm{V}\cdot\nabla\phi)-gW]=0$ (5)

Scalar moisture equations

$\partial_{t}Q_{m}+(\nabla\cdot\bm{V}q_{m})=F_{Q_{m}}$ (6)

where

$\nabla\cdot\bm{V}a=\partial_{x}(Ua)+\partial_{y}(Va)+\partial_{\eta}(\Omega
a)$

and

$\bm{V}\cdot\nabla a=U\partial_{x}a+V\partial_{y}a+\Omega\partial_{\eta}a,$

for a generic variable $a$.

These equations are coupled with the diagnostic equation for dry hydrostatic
pressure

$\partial_{\eta}\phi=-\alpha_{d}\mu_{d}$ (7)

and the  diagnostic relation for the full pressure (dry air combined with
water vapor)

$p=p_{0}\Big{(}\frac{R_{d}\theta_{m}}{p_{0}\alpha_{d}}\Big{)}^{\gamma}$ (8)

In Eq. (2) - (6), $\alpha_{d}=\frac{1}{p_{d}}$ is the inverse density of dry
air, $\alpha=\alpha_{d}(1+q_{v}+q_{c}+q_{r}+...)^{-1}$, while $R_{d}$ is the
gas constant for dry hair and $\gamma=\frac{c_{p}}{c_{v}}=1.4$ is the ratio of
the heat capacities for dry air. The right hand side terms
$F_{U},F_{V},F_{W},F_{\Theta_{m}}$ contain the Coriolis and curvature terms
along with mixing terms and physical forcings.

Since the ARW solver works on a grid, a projection from the computational flat
space to the physical spherical space is needed to interpret results: in
particular, ARW implements the projection using map scale factors $m_{x}$ and
$m_{y}$, defined as the ratio of the distance ($\Delta x,\Delta y)$ in
computational space to the corresponding distance on the earth’s surface:

$(m_{x},m_{y})=\frac{(\Delta x,\Delta y)}{\text{distance on the earth}}.$

In practice, the ARW actually solves a perturbative form ([17], pp.12) of the
previous governing equations, where the momentum variables are redefined as
follows

$U=\frac{\mu_{d}u}{m_{y}},V=\frac{\mu_{d}v}{m_{x}},W=\frac{\mu_{d}w}{m_{y}},\Omega=\frac{\mu_{d}\omega}{m_{y}}.\\\
$

To solve the governing equations, ARW adopts a time-split integration scheme.
In atmospheric model, however, slow and fast processes coexist and they have
to be integrated differently: the former through a third-order Runge Kutta
(RK3) time integration scheme; while the latter over a smaller RK3 time steps
to guarantee numerical stability.

#### Level set equation

In SFIRE [18], the propagation of a fire burning in the area
$\Omega=\Omega(t)$ in the horizontal $(x,y)$ plane on which the Earth is
projected, is implemented by the level set method, which evolves a function
$\psi=\psi(x,t)$, called the level set function, such that the burning area at
a time $t$ is

$\Omega(t)=\\{x:\psi(x,t)\leq 0\\}$

and the fireline $\Gamma(t)$, i.e. the boundary of the burning region
$\Omega(t)$, is the level set

$\Gamma(t)=\\{x:\psi(x,t)=0\\}.$ (9)

The model adopts a semi-empirical approach, where the fire spread rate $S$ is
computed from fuel properties, using the following modified Rothermel formula
[19]:

$S=R_{0}(1+\phi_{W}+\phi_{S}).$ (10)

Here, $R_{0}$ represents the spread rate in the absence of wind, whereas
$\phi_{W}$ and $\phi_{S}$ are respectively the wind factor and the slope
factor. Given a point $x\in\Omega(t)$, the time of ignition $t_{i}$ is defined
as the time when the point $x\in\Gamma(t_{i}(x))$. The model assumes that at
each location on the grid the fuel fraction, at the beginning $F=1$, decreases
exponentially from the time of ignition $t_{i}$.

On the fireline, from Eq. (9), the tangential component of the gradient
$\nabla\psi$ is zero, meaning that the fire propagation speed is normal to the
fireline $\Gamma$: the model postulates $S$ is a function of the normal
component of the wind factor $U$ and the terrain gradient $\nabla z$ [19].

The evolution of the level set function is governed by the partial
differential equation

$\partial_{t}\psi+S||\nabla\psi||=0,$ (11)

called the level set equation, which is solved numerically through
discretization on the refined fire grid, adopting Heun’s method, a second-
order Runge Kutta method ([18], p.596).

###  3.6 The new architecture

Figure 3: Physics-informed neural network architecture for solving PDEs

We have replaced the original numerical solvers with a profoundly different
approach, by taking advantage of the recently developed Julia packages
NeuralPDE.jl and DiffEqFlux.jl, which support Physics Informed Neural Networks
(PINNs) for automated PDE solving and Backwards Stochastic Differential
Equation (BSDE) methods to deal with parabolic PDEs. We employed the former
method: PINNs are feed forward neural networks that are trained to solve
supervised learning tasks while respecting any given law of physics described
by general nonlinear partial differential equations. The resulting neural
networks form a new class of universal function approximators that naturally
encode any underlying physical laws as prior information [20].

The underlying architecture can be summarized by the following steps, with
reference to Fig.3. Let’s put ourselves in the most general case, where we
want to find n target functions $[u_{1},...,u_{n}]$ which satisfy a system of
n PDEs. We first construct a surrogate solution to u(x) as a neural network
û(x;$\theta$) with n inputs and outputs and parameters $\theta$; since a
neural network is mathematically a composite function, the derivatives of û
with respect to its input can be evaluated by applying the chain rule for
differentiating compositions of functions using automatic differentiation
(AD). Then we need to define a loss function (we used the $L^{2}$ norm) in
order to calculate the discrepancy between $[u_{1},...,u_{n}]$ and the outputs
of the neural network: this is where we force our network to satisfy the
physics imposed by the PDEs. Indeed, the loss is defined by substituting our
neural network and its derivatives back into the equations we want to solve
(where we brought all the terms on the left-hand side of the equal sign). The
last step is to train the neural network to find the best parameters by
minimizing the loss function with gradient-based optimizers, such as Adam or
LBFGS. It is straightforward that minimizing the loss function means solving
the equations, since we are approximating our target functions with û better
and better. The key point is that it converts an integration problem into a
mere minimization task and it does not need data to train, because it is
trained on the PDEs themselves.

###  3.7 Model selection process that led us to this choice

Our aim from the beginning has been finding an ML alternative to standard
numerical methods. However, we have explored different architectures within
the field of Scientific Machine Learning. We started our investigations using
the DiffEqFlux.jl library [21], which defines and solves neural ordinary
differential equations (i.e. ODEs where a neural network defines its
derivative function), which are implemented in the following way. First we
need to get data from the numerical ODE solver, then define a neural network
with a NeuralODE layer, which constructs a continuous-time recurrent neural
network: at a high level this corresponds to solving the differential equation
during the forward pass and using a second differential equation that
propagates the derivatives of the loss function backwards in time. Now it’s
possible to train the model using an optimizer (ADAM or LBFGS) which minimizes
the loss function. We have consulted the library documentation and the GitHub
resources to implement some toy models with this approach, such as the Lotka-
Volterra equation, the 2D reaction diffusion equation, the 1D Fisher-KPP and
the linear Burgers equation. Furthermore we tried to go beyond these
implementations: we solved the reaction-diffusion equation with a
Convolutional Neural Network (CNN) to investigate peculiarities of a CNN with
respect to a common feed forward network in approximations of diffusive
models. When solving the Fisher-KPP we ran into a problem: since the neural
ODE approximates the derivative of the solution, after it is trained we have
not reached the desired solution yet; we still have to use a numerical solver.
Following our initial purpose, we chose to add a Recurrent Neural Network on
top of the NN + CNN architecture used in the first step. Clearly this model
showed some relevant disadvantages:

  * •

The predicted slices were very noisy, because the integrator would have
required a greater amount of data for training;

  * •

The integrator tended to overfit during the training, getting a model not
applicable on other time scales and too problem-specific.

This kind of issue made us lean towards the NeuralPDE PINNs approach, that is
our ultimate architecture, clearly described in paragraph IV-D where the NN
directly approximates the objective function, without needing further
integrations. Before attempting to model the complex PDEs inside WRF, we
characterized the quality of the architecture by running it on the linearized
Burger equation, non linear Burger equation without viscosity, non linear
Burger equation in the general form and the Poisson equation in 2 dimensions.

###  3.8 Implementation of the Level set

The level-set differential equation is, as previously explained, the
mathematical core for calculating the spread of the fire. The level set
equation (that is a 2D surface embedded in a 3D euclidean space) is initially
set as the distance of each point from the ignition’s line. In our approach it
is required that every quantity is provided either as a constant or as a
function. Therefore we implemented the initial condition as a cone with
elliptical section, because we had to deal with circular and linear ignition
shapes (theoretically, any ignition shape can be given in input). The cone is
then shifted down in order to make the contour level $z=0$ resembling the
desired initial fire line. Multiple ignition points are possible, but further
refinements are needed to make this possibility fully functional. The
algorithm contains all the necessary expressions for the calculation of the
fire spread rate, taking into account the values of the wind, the altimetric
profile, the fuel map etc. These variables should be given in input as
differentiable vector fields of space and time. This can be considered both a
limit and an advantage, in fact on one hand the discretized measurements must
be fitted using approximating functions, but on the other hand this method is
very computationally efficient. The possibility to give in input finite matrix
of values is presented in "Future Work" and depends on the fact that some key
libraries we use are still under development. You can see in the appendix an
example of the altimetric profile of the Isom Creek fire location fitted using
a 4th degree polynomial. The last step before the training of the model is a
proportional scaling of the values to a domain with size in the order of
$10^{1}$ and the definition of a discretization that is used only to run the
training algorithm. The outputs will be continuous functions. The
discretization to this interval prevents some numerical instabilities and
spikes that are linked to the early stage of development of some Julia
libraries in use. The loss function is made of two parts: the former
minimizing the $L^{2}$ error referred to the PDE definition and the latter
minimizing the $L^{2}$ distance between the boundary conditions and the
solution evaluated at the boundaries. At this point the training algorithm is
instantiated and launched. The minimization of the loss functions is the
process that actually solves the PDE and constitutes the main load for the
CPU. It can be easily parallelized and thus accelerated using GPUs. When the
training is completed the prediction undergoes a new proportional scaling that
gives back the original domain shape.

###  3.9 Implementation of the Euler system

Writing and solving the 7-equation Euler system in Julia was really
challenging, in fact at the moment we are not aware of any publications where
these techniques are yet applied to PDE systems of such complexity. The first
problem was to figure out which were the independent and dependent variables.:
it might seem rather simple, but climate models adopt different conventions
depending on the circumstances (e.g., they generally use pressure coordinates
instead of height coordinates). That is why we presented various Julia files,
both with six and seven target functions, for reasons we cannot explain here
for brevity. Since the chain rule for derivatives hasn’t been implemented into
ModelingToolkit.jl yet, one of the most demanding operations was writing all
the derivatives involved by symbolic calculation. Then we had to choose the
most suitable training strategy, which turned out to be QuadratureTraining. We
have used parabolic initial conditions for a first evaluation, but we want to
investigate them more in depth, since the study of boundary conditions often
deserves a paper on its own. Unfortunately, the NeuralPDE library is still
unable to treat this kind of problem with stability, and often incurs errors
due to internal divergence calculations. Despite this, we have been able to
obtain convergence of the loss function, although it is not enough to present
valid results. We contacted the authors of these libraries, that are still
under development, and we are looking forward to contributing.

###  3.10 Optimization and coding

For the implementation of the model, we chose to rely on the Julia programming
language. Julia, despite being a very young language with a huge amount of
bugs and improvements still to be implemented, compared to the more diffused
Python, offers better performance since it is not an interpreted language
providing Just in Time (JIT) compilation with different optimization levels
[1]. Moreover, Julia is particularly well suited for numerical computation and
for the solution of complex physical models, both considering the syntax and
the presence of some specific libraries such as DifferentialEquations.jl and
ModellingToolkit.jl.

The model was implemented using the low-level interface of the NeuralPDE.jl
library which contains the necessary methods for the generation of the
training datasets and of the loss functions starting from the explicit form of
the equations and the boundary conditions. This library allows some very high
level macros for the declaration of derivatives, functions and parameters and
allows you to choose the differentiation techniques (numerical or automatic)
and the optimization engine to be used at a later time. The library relies in
turn on GalacticOptim.jl, which provides the methods for minimizing the loss,
and on ModelingToolkit.jl, which provides the constructors for the derivatives
and on Flux.jl, which is the library that implements the constructors of
neural networks with the respective backpropagation and prediction functions.

###  3.11 Synthetic datasets - ideal simulations

The WRF source code includes some idealized simulations that are used as a
playground for the model: these simulations are composed of an namelist.input
file and an input_sounding file. The former contains both data about the
physics models that should be used, as well as terrain information, to
generate the surface height mesh; the other is a standard way of storing
weather information, containing potential temperature, water vapor mixing
ratio, and U and V wind speed components values at various pressure values.

We decided that it was best to start with an ideal simulation to be able to
control the response of the model, so we modified the already present
two_fires test to fit our needs, reducing the number of fires to one and
increasing the simulation time to 60 minutes.

The terrain generated is completely flat, and the entire domain is covered in
tall grass (Anderson fuel category 3 [22]). Since the wind is entirely along
the V component, we expect the fire to expand upwards.

###  3.12 Real datasets - real simulations

In order to run a real data wildfire simulation in the WRF model, static
terrestrial data, atmospheric data and a map of the 13 Anderson’s fuel
category of the desires domain are required. The geographic global static data
containing the mandatory fields requested are provided by the National Center
for Atmospheric Research (NCAR) website [23]; whereas an Anderson 13 fuel
categories [22] map and a high-resolution topographic data (about 5 meters
resolution) are provided by the Landfire [24] site. Lastly, for our simulation
we adopted historical meteorological fields from a Global Forecast System
(GFS) reanalysis, with a resolution of 0.25 degree and provided with a 6-hour
time step. The latter are accessible in Grib2 format on the Research Data
Archive RDA of NCAR website [25].

All the datasets were appropriately pre-processed by the WRF Pre-Processing
System (WPS), which allows to define the simulation domain. In addition,
through a horizontal interpolation of datasets, WPS provides the initial and
boundary conditions that will be used by the WRF model. For this work we
carried out a simulation of the Isom Creek fire, Alaska (USA), [26] where the
ignition happened on 05 June 2020 at 03:15 PM and the progression lasted until
13 June, reaching an extension area of 12,180 acres (Fig.4).

Figure 4: Isom Creek Fire progression map Figure 5: Altimetric profile
interpolation of the Isom Creek fire region.

Therefore, the center of the simulation domain has been fixed at the
coordinates of the Isom Creek fire ignition point (65.854; -149.579). The
chosen atmospheric domain is a mesh of 97x97 grids, each with a side of 100 m,
while the 2D mesh of the fire domain is 20 times finer, therefore with a
resolution of 5 meters. The extrapolation of the domain from the land-surface
datasets takes place by means of the geogrid.exe program in WPS, and at this
point occurs also the appropriate overlap between all the input geographical
data Fig.6 and Fig.7

Figure 6: Comparison between global topological data and high resolution data
after geogrid.exe program run Figure 7: Overlap of the 13 Anderson’s
categories map on the domain of simulation after the geogrid.exe program run

By briefly analyzing these datasets we found a few interesting facts. The most
noticeable aspect is the fuel distribution: the most abundant fuel is timber
(category 10), followed by brush (category 5) and closed timber litter
(category 8). This fuel composition is wide-spread in this area and due to the
presence of timber it can lead to potential fire control difficulties [22].
The presence of a river in the top part of the domain suggests that the fire
will not spread here.

The mountainous topography of the area will also affect the spread of the
fire, as the terrain gradient is part of the calculation for the fire spread
rate.

ungrib.exe deals with the pre-processing of global atmospheric dynamical data.
Finally, the metgrid.exe program is responsible of the combination of
meteorological data to the simulation domain.

Some attempts of simulations with a second finer domain, inner to the domain
described above, have been made. The side of the inner domain is a third of
the outer domain, and the number of cells is also 97, therefore we obtained a
domain of about 3 km with a resolution three times higher. From these tests we
obtained relatively valuable simulations, bearing in mind that for larger
simulations researchers set the beginning of the simulation at least 12 hours
before of the period of interest. However, we preferred to come back to a
single domain simulation in order to speed up the run, paying the price in
terms of a meteorological simulation affected by boundary effects. Our
ultimate Isom Creek fire simulation ran for 18 hours, starting from the 5th
June at 12:00 PM.

###  3.13 Ideal simulation: "One fire"

For clarity, we renamed the Two Fires ideal simulation as One Fire, after
having removed one ignition point. It is a simple case that demonstrates that
our model works well in case of fires that spread in flat areas where the fuel
doesn’t change. It also evolves according to the wind direction, represented
by a vector with two components. The results at different time frames are
represented in Fig.8. The details about the architecture employed and the
training are listed in TABLE 2.

Optimization Algorithm |  ADAM
---|---
Iterations |  4800
Final Objective Value |  6.39 e-8
Training Strategy |  QuadratureTraining()
Domains |  $t\in[0,10]$, $x\in[0,10]$, $y\in[0,10]$
Training Mesh Size |  $[dt,dx,dy]=[0.17,0.02,0.02]$
Boundary Condition |  $u(0,x,y,\theta)=((5(x-0.3))^{2}+(0.15y)^{2})^{\frac{1}{2}}-0.2$
Neural Network dimensions |  $3>16>1$
Training Time |  647 s
Table 2: Technical Specs of the model - One Fire Figure 8: Comparison between
the WRF and PINNs outputs at different time frames

In order to provide a quantitative measure of the error between the outputs we
decided to use the Hausdorff distance [27],

$d_{\mathrm{H}}(X,Y)=\max\left\\{\,\sup_{x\in X}\inf_{y\in
Y}d(x,y),\,\sup_{y\in Y}\inf_{x\in X}d(x,y)\,\right\\},$

where $X$ and $Y$ are the firelines to be compared. This distance is further
normalized on the area of the curves that represent the fireline, as Fig.9
shows. The fact that the error decreases tells us that as the area grows with
time, the difference between the two regions doesn’t increase, meaning that
they evolve in synchrony and with the same shape.

Figure 9: Hausdorff distance as a measure of the error - One Fire.

###  3.14 Simulation of a real wildfire: Isom Creek, Alaska, 6/5/2020

We have been able to simulate the Isom Creek fire both with WRF and with our
implementation. Although it is not possible to perform a precise quantitative
comparison between WRF and the real data, or between our result and the real
data (due to a lack of temporal precision of the available data at [28]), it
is clear that the WRF output simulates the wildfire accurately, at least
visually (see Fig.4). The result obtained with the PINN architecture is shown
in Fig.10. It is clear that our implementation hasn’t been able to capture the
shape of the fireline with the same accuracy. This will be pointed out below
as a temporary limitation to our model: at the moment it is not possible to
assign a different fuel category to each mesh point (i.e. using a matrix to
build a function), so the fireline will always have a symmetrical and smooth
perimeter. The specifics of the architecture follow.

Optimization Algorithm |  ADAM
---|---
Iterations |  3000
Final Objective Value |  2.14 e-6
Training Strategy |  QuadratureTraining()
Domains |  $t\in[0,10]$, $x\in[0,10]$, $y\in[0,10]$
Training Mesh Size |  $[dt,dx,dy]=[0.17,0.02,0.02]$
Boundary Condition |  $u(0,x,y)=((x-5)^{2}+(0.7(y-5))^{2})^{\frac{1}{2}}-0.2$
Neural Network dimensions |  $3>16>1$
Training Time |  270 s
Table 3: Technical Specs of the model Figure 10: Comparison between the WRF
and PINNs outputs at different time frames

The error can be quantitatively measured with the Hausdorff distance as
before, with the outcome plotted in Fig.11.

Figure 11: Hausdorff distance as a measure of the error - Isom Creek.

The considerations outlined in the previous paragraph remain valid here (see
10 for additional comparison plots).

##  4 DISCUSSION

###  4.1 Discussion of the experimental results

The results exposed above show the feasibility of utilizing the Julia
NeuralPDE library to model a wildfire related differential equation, i.e. the
level set equation. As disclosed above, we have also obtained some preliminary
but unstable results regarding the atmospheric coupling of the model by
attempting to solve the Flux-Form Euler system equation. The outputs obtained
with our approach don’t strictly follow the jagged firelines of WRF: that’s
not an intrinsic limit of the architecture, which doesn’t lack variance even
if it is rather simple; the reason is, as explained above, that the inputs
(fuel and wind in particular) we gave to the NN are not varied enough, due to
limitations of the library. In spite of this we have shown the ability of our
implementation to reproduce the desired time evolution of the fireline as
given by WRF.

###  4.2 Advantages of this architecture

The various advantages of this approach are:

  * •

We can get a rough estimate of the speed up our model provides - compared to
numerical solvers - by noting that the train time is of some hundreds of
seconds, while the WRF run time is more than an hour. This point has to be
further studied: on one hand, we used a simpler set of inputs for the reasons
we will explain later; on the other hand we worked on one of the eight cores
at our disposal.

  * •

While numerical methods only solve on the initial domain range, this
architecture learns the neural network parameters that can be used to extend
the predicted solution outside of the training domains (with a decrease in
accuracy if the extrapolation is extended too much). This can be used for
"forensic investigations": it is possible to reconstruct the ignition point
and the evolution of the fireline from the final state of the wildfire.

  * •

It is possible to interpolate the solution on a continuous mesh.

  * •

Modifying the equations of the model is as easy as changing a few lines of
code, instead of re-implementing the discretization and perturbative form of
the equations, allowing for faster investigation of new models.

  * •

Given the possibility to use several CPU cores, the speed could increase even
more and this approach could be run iteratively in order to simulate the
outcome of different containment scenarios and choose the best among them.

  * •

This architecture doesn’t need neural networks with a high number of degrees
of freedom to be accurate, reducing the overall computational cost and
increasing performance. In quantitative terms, this means not exceeding some
tens of neurons on a single hidden layer.

  * •

As a consequence, the model is more interpretable (it is less similar to a
complete black-box) and its users can run the resulting simulator on standard
PCs without needing high performance machines.

  * •

PINNs can be used to extrapolate physical laws once they have been
interpolated, hence providing a tool to refine theoretical models. However,
this requires a greater amount of input data.

###  4.3 Limitations

Let’s now analyze the limits of our approach:

  * •

At the moment, it cannot be run on GPUs nor multiple CPUs because the library
we used is not fully developed yet. We contacted the developers and they are
working on it.

  * •

The library doesn’t support the usage of matrices as functions of their
indices. This means it is not possible to use a variable fuel map, fundamental
in order to capture the evolution of the fire with higher accuracy. We opened
an issue about this [29].

  * •

Our approach suffers from numerical instability depending on the training
strategy adopted and on the number of neurons and layers.

  * •

The model only works on relatively small domains and is unstable for larger
ones. Therefore we had to scale down all the quantities involved and introduce
scale factors on the fire spread rate in order to make it compatible to the
evolution of the output of WRF.

  * •

We still have to figure out how to predict the spread of a wildfire with more
than one ignition point.

  * •

The architecture has not yet been extended to Convolutional Neural Networks,
Recurrent Neural Networks or others.

##  5 FUTURE WORK

This research is the first step towards a concrete implementation of PINNs-
based wildfire simulation. We have designed a chart (Fig.12) that shows the
next improvements to be made in order to achieve such a ambitious goal.

Figure 12: Future prospect of our work

##  6 HARDWARE SUPPORT AND CODE

The implementation described in the previous paragraph is based on the ML
framework Flux.jl that is agnostic with respect to the array type. This allows
to also pass to this engine GPU arrays (CUDA arrays) that are allocated and
computed on GPUs (General Purpose GPUs). For the development we used a
workstation with the following technical specifications:

  * •

CPU: Intel Core i7 10700K 8 cores / 16 threads

  * •

RAM: 64 GB DDR4 unbuffered

  * •

GPU: Nvidia 2080 Ti 11 GB GDDR6

Unfortunately, the GPU interface of NeuralPDE.jl is not mature yet.

All the code produced is published on a public repository available at this
link:
https://github.com/MachineLearningJournalClub/MLJC-UniTo-ProjectX-2020-public.
The output of the WRF simulations can be found in the following Google Drive
folder:
https://drive.google.com/drive/folders/1wUCKUyVwC0Pf-e9WlLiqOxRLF0or2D0U?usp=sharing
or https://tinyurl.com/mljc-unito-px2020.

##  7 SOCIETAL CONSEQUENCES

###  7.1 The impact of real time simulations for containment

As previously mentioned, one of the most important aspects of implementing a
fast simulator is the possibility of employing it in wildfire containment, as
firefighters could simulate the possible outcomes for different strategies,
therefore applying the best one. At the moment WRF-SFIRE is not fast enough to
be used this way, but the ML approach shows promising results for such an
application. As of now, firefighters still struggle to predict changes in fire
spread, as we observed even in the Isom Creek fire reports [30], which lead to
the destruction of over 12,000 acres of forest, a huge economic cost, and the
damage done to an entire ecosystem. The usage of simulators with the
characteristics we outlined would have the ethical value of saving lives, on
top of saving huge amounts of money by reducing the targeted area.

###  7.2 Simulation of fire propagation and climate change

A fast and reliable wildfire simulation can also be used to simulate how
wildfires might develop and worsen in different climatic conditions, for
example by reproducing the same fire in various scenarios. This can be done by
changing some of the input meteorological data, such as the temperature and
the radiative forcing (RF). Using predictions of how climate might change due
to global warming, the model could prove the impact that this crisis might
have on wildfires.

###  7.3 Forensic investigation

A consequence of how boundary conditions are handled in our model is that,
provided the fire shape in a certain instant, it would be possible to
extrapolate the fireline in both "directions" of time. There is usually very
little information about the temporal evolution of the fire (especially in
remote areas), so the ability to extrapolate the shape of the spread backwards
in time is crucial to understand where and how the fire started. These
"forensic" investigations can provide insights that are valuable for fire
prevention and protection. Our model is able to do so without any addition.
The only requirement is that a "reasonable" level set is defined: as for the
initial value problem, the level set equation must be set such that it does
not have too big or too small values. We suggest to keep it in the range
$(-20,20)$, however this is not necessary and is situation-dependent. It is
very frequent that the negative part of the level set equation is indeed in
the range $(-1,0)$. An example is portrayed in Fig.13.

Figure 13: Example of backward time prediction

##  8 CONCLUSIONS

The study we carried out had the goal to investigate the applicability of the
recently developed field of Scientific Machine Learning on climate, wildfire
in particular, models. We have outlined some results that tell us that many
improvements are needed in order to transform this into a validated product,
but also show the big potential of our approach. We need to add further
refinements to the implementation in order to carry out a precise time
comparison between the two approaches, but the results obtained thus far show
promising evidence. The encouraging outcome inspires us to continue our work
by improving the architectures and possibly employ them in different fields of
research. We hope that this line of research will be considered as a starting
point for a more effective cohesiveness between Machine Learning and Physical
Models in Climate Science, and thus further explored by other researchers.

##  9 ACKNOWLEDGEMENTS

This work was presented at the ProjectX 2020 competition by UofT AI. We
acknowledge University of Turin, Machine Learning Journal Club for supporting
us. We thank Professor Enrico Ferrero (Università del Piemonte Orientale),
Professor Massimiliano Manfrin (University of Turin) and the whole Atmospheric
Physics and Metereology Group, PhD Christopher Rackauckas (Massachussets
Institute of Technology), PhD Kirill Zubov (Saint-Petesburg State University),
Vaibhav Dixit (Julia Computing), PhD Brian Wee (Founder at Massive
Connections), Dr. Rustem Arif Albayrak (NASA), PhD David Marvin (CEO at Salo
Sciences), Professor Piero Fariselli (University of Turin) and Pietro
Monticone M.Sc. student (University of Turin), for their precious help and
availability. We acknowledge the company Mollificio Astigiano (Belveglio,
Asti, Italy) for providing the computational power needed for this research
and the HPC4AI center of the University of Turin for their support.

##  10 APPENDIX

Figure 14: Perimeter-normalized Hausdorff distance as a measure of the error
between outputs - One Fire Figure 15: Perimeter-normalized Hausdorff distance
as a measure of the error between outputs - Isom Creek Figure 16: Time
evolution of the perimeters of the firelines - One Fire Figure 17: Time
evolution of the perimeters of the firelines - Isom Creek Figure 18: Time
evolution of the areas of the firelines - One Fire Figure 19: Time evolution
of the areas of the firelines - Isom Creek Figure 20: Boundary condition
comparison and training - One Fire Figure 21: Boundary condition comparison
and training - Isom Creek

## References

  * [1] Miles Lubin and Iain Dunning  “Computing in Operations Research Using Julia”  In _INFORMS Journal on Computing_ 27.2 Institute for Operations Researchthe Management Sciences (INFORMS), 2015, pp. 238–248  DOI: 10.1287/ijoc.2014.0623
  * [2] Topher Gauk-Roger, Stella Chan, Jason Hanna and Steve Almasy  “California wildfires: Fire chief says dozens of major blazes have state in ’dire situation”’  In _CNN_ , 2020
  * [3] “Fire statistics, CALFIRE”  URL: https://www.fire.ca.gov/incidents/2020/
  * [4] Marco Turco, Maria-Carmen Llasat, Jost Hardenberg and Antonello Provenzale  “Climate change impacts on wildfires in a Mediterranean environment”  In _Climatic Change_ 125.3-4 Springer ScienceBusiness Media LLC, 2014, pp. 369–380  DOI: 10.1007/s10584-014-1183-3
  * [5] David M. Romps, Jacob T. Seeley, David Vollaro and John Molinari  “Projected increase in lightning strikes in the United States due to global warming”  In _Science_ 346.6211 American Association for the Advancement of Science, 2014, pp. 851–854  DOI: 10.1126/science.1259100
  * [6] Alejandra Borunda  “The science connecting wildfires to climate change”  URL: https://www.nationalgeographic.com/science/2020/09/climate-change-increases-risk-fires-western-us/
  * [7] “Impacts, Risks, and Adaptation in the United States: The Fourth National Climate Assessment, Volume II”  U.S. Global Change Research Program, 2018  DOI: 10.7930/nca4.2018
  * [8] Bruce Lieberman  “Wildfires and climate change: What’s the connection?”  URL: https://yaleclimateconnections.org/2019/07/wildfires-and-climate-change-whats-the-connection
  * [9] Christopher Rackauckas et al.  “Universal differential equations for scientific machine learning”  In _arXiv preprint arXiv:2001.04385_ , 2020
  * [10] Věra Kůrková  “Kolmogorov’s theorem and multilayer neural networks”  In _Neural Networks_ 5 Elsevier Science, 1992  DOI: 10.1016/0893-6080(92)90012-8
  * [11] Yu Li et al.  “Deep learning in bioinformatics: Introduction, application, and perspective in the big data era”  In _Methods_ 166 Elsevier, 2019, pp. 4–21
  * [12] Luca Cenci, Luca Pulvirenti, Giorgio Boni and Nazzareno Pierdicca  “Defining a Trade-off Between Spatial and Temporal Resolution of a Geosynchronous SAR Mission for Soil Moisture Monitoring”  In _Remote Sensing_ 10.12 MDPI AG, 2018, pp. 1950  DOI: 10.3390/rs10121950
  * [13] Maziar Raissi, Paris Perdikaris and George E Karniadakis  “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations”  In _Journal of Computational Physics_ 378 Elsevier, 2019, pp. 686–707
  * [14] Philipp Grohs, Arnulf Jentzen and Diyora Salimova  “Deep neural network approximations for Monte Carlo algorithms”, 2019  arXiv:1908.10828 [math.NA]
  * [15] T. Gleixner and I. Molnar  “perf Wiki Homepage”, https://perf.wiki.kernel.org/index.php/Main_Page
  * [16] “User’s Guides for the Advanced Research WRF (ARW) Modeling System, Version 4”  URL: https://www2.mmm.ucar.edu/wrf/users/docs/user_guide_v4/contents.html
  * [17] William C. Skamarock et al.  “A Description of the Advanced Research WRF Model Version 4”  UCAR/NCAR, 2019  DOI: 10.5065/1DFH-6P97
  * [18] J. Mandel, J.. Beezley and A.. Kochanski  “Coupled atmosphere-wildland fire modeling with WRF 3.3 and SFIRE 2011”  In _Geoscientific Model Development_ 4.3, 2011, pp. 591–610  DOI: 10.5194/gmd-4-591-2011
  * [19] Richard C. Rothermel  “A Mathematical Model for Predicting Fire Spread in Wildland Fires”  In _USDA Forest Service Research Paper INT-115_ , 1972  URL: http://www.treesearch.fs.fed.us/pubs/%2032533
  * [20] Kirill Zubov  “Physics-informed neural networks (PINNs) solver on Julia. GSoC 2020. First evaluation.”, 2020  NEXTJOURNAL: https://nextjournal.com/kirill_zubov/physics-informed-neural-networks-pinns-solvers-on-julia-gsoc-2020-first-evaluations
  * [21] Christopher Rackauckas et al.  “DiffEqFlux.jl - A Julia Library for Neural Differential Equations”  In _CoRR_ abs/1902.02376, 2019  arXiv: http://arxiv.org/abs/1902.02376
  * [22] Hal E. Anderson  “Aids to determining fuel models for estimating fire behavior”  U.S. Department of Agriculture, Forest Service, Intermountain ForestRange Experiment Station, 1982  DOI: 10.2737/int-gtr-122
  * [23] “WPS V4 Geographical Static Data”  URL: https://www2.mmm.ucar.edu/wrf/users/download/get_sources_wps_geog.html
  * [24] “Landscape Fire and Resource Management Planning Tools (LANDFIRE)”  URL: https://landfire.gov/
  * [25] “NCEP GFS 0.25 Degree Global Forecast Grids Historical Archive”  URL: https://rda.ucar.edu/datasets/ds084.1/
  * [26] “Isom Creek Fire Information”  URL: https://inciweb.nwcg.gov/incident/6758
  * [27] Christian Knauer, Maarten Löffler, Marc Scherfenberg and Thomas Wolle  “The directed Hausdorff distance between imprecise point sets”, 2009  arXiv:0909.4642 [cs.CG]
  * [28] “Isom Creek Fire Progression Map”  URL: https://inciweb.nwcg.gov/incident/map/6758/7/98095
  * [29] “Embedding arrays into NeuralPDE differential equations”  URL: https://github.com/SciML/NeuralPDE.jl/issues/177
  * [30] “AK Fire Info - Isom Creek Fire Reports”  URL: https://akfireinfo.com/tag/isom-creek-fire/