# Using Echo State Networks to Inform Physical Models for Fire Front
Propagation

Myungsoo Yoo
University of Missouri
<EMAIL_ADDRESS>
&Christopher K. Wikle
University of Missouri
<EMAIL_ADDRESS>
Author of correspondence

###### Abstract

Wildfires can be devastating, causing significant damage to property,
ecosystem disruption, and loss of life. Forecasting the evolution of wildfire
boundaries is essential to real-time wildfire management. To this end,
substantial attention in the wildifre literature has focused on the level set
method, which effectively represents complicated boundaries and their change
over time. Nevertheless, most of these approaches rely on a heavily-
parameterized formulas for spread and fail to account for the uncertainty in
the forecast. The rapid evolution of large wildfires and inhomogeneous
environmental conditions across the domain of interest (e.g., varying land
cover, fire-induced winds) give rise to a need for a model that enables
efficient data-driven learning of fire spread and allows uncertainty
quantification. Here, we present a novel hybrid model that nests an echo state
network to learn nonlinear spatio-temporal evolving velocities (speed in the
normal direction) within a physically-based level set model framework. This
model is computationally efficient and includes calibrated uncertainty
quantification. We show the forecasting performance of our model with
simulations and two real data sets - the Haybress and Thomas megafires that
started in California (USA) in 2017.

_K_ eywords wildfire $\cdot$ signed distance function $\cdot$ boundary $\cdot$
level set method $\cdot$ echo state network $\cdot$ spatio-temporal

##  1 Introduction

Wildfires can negatively affect natural and human ecosystems through loss of
life, property damage, emission of air pollutants and impacts on water
quality, among others (e.g., Brotons et al.,, 2013; Hohner et al.,, 2019).
Indeed, wildfires can continue to pose a risk even after they are
extinguished. For instance, Cannon and DeGraff, (2009) raised public awareness
of the impact of the post-fire debris flow. As stake holders consider
tradeoffs between managing ecosystems and preventing loss of life and
property, having as much information to plan wildfire suppression tactics as
early as possible is important. This requires an efficient and reliable
forecasting model for the evolution of wildfires.

The level set method proposed by Osher and Sethian, (1988) has become one of
the most useful approaches for modeling the evolution of wildfire fronts in
operational fire models (e.g., Miller et al.,, 2015). Osher and Fedkiw, (2002)
summarized various implementations of the level set method, and Gibou et al.,
(2018) provided a recent review. In the application to wildfire, the level set
method frequently relies on heavily-parameterized quasi-empirical formulations
of the wildfire front velocity, such as the Rothermel model (Rothermel,,
1972). Such models can be suitable and perform well under homogeneous
environmental and atmospheric conditions. For example, Dabrowski et al., 2022a
used the Rothermel model with an unknown parameter to inform velocity in the
level set method by using an ensemble Kalman filter (EnKF). Similarly, Mallet
et al., (2009) and Alessandri et al., (2021) simplified the velocity through
formulas proposed by Fendell and Wolff, (2001) and Lo, (2012), respectively.
None of these approaches utilizes a data-driven approach to learn the
spatially and temporally varying front velocity over the entire domain, which
is required for realistic front propagation in complex environments,
especially where large fires create their own “weather”.

As discussed below, the implicit function used in the level set method to
represent the fire boundary is a time-varying signed distance function that is
referenced in space and time. Spatio-temporal data are ubiquitous in the
environmental sciences and multi-level hierarchical representations have long
been used to model such complex processes (e.g., see Cressie and Wikle,, 2011;
Banerjee et al.,, 2015, for general discussion). Likewise, the dependency
observed in space or time has triggered intellectual curiosity in the deep
learning community. Wikle and Zammit-Mangion, (2022) provide a recent overview
of statistical and AI approaches for deep learning spatial and spatio-temporal
data. Convolutional neural networks (CNNs), recurrent neural networks (RNNs),
and their variants are the primary approaches in the deep learning literature
for such data. In particular, hybrid convolutional recurrent networks (e.g.,
Zuo et al.,, 2015) have worked well for spatio-temporal data. For example, in
the context of wildfire, Burge et al., (2022) proposed the EPD-convLSTM model,
which accounts for spatial and temporal dependency by combining a
convolutional long short-term memory (LSTM) RNN and an Encoder-Processor-
Decoder (EPD). However, these methods do not provide a mechanism for
uncertainty quantification, do not directly incorporate physical relationships
(such as implied by the level set method), and require a large amount of
training data.

An enhancement to traditional deep learning neural network approaches is to
use so-called “physics-informed neural networks” that seek to encode physical
relationships into the deep neural architecture (e.g., Raissi et al.,, 2019).
These methods have recently been considered for level-set fire front
propagation in Dabrowski et al., 2022b . Their work is implemented in a
Bayesian framework that allows for data assimilation and uncertainty
quantification. However, this approach does not consider realistic data-driven
spatial and temporally varying velocities or large-scale megafire spread,
which is our interest here. Alternatively, Yoo and Wikle, (2022) considered a
hybrid physical/statistical level set model within a hierarchical Bayesian
framework. Their approach included a latent vector autoregressive stochastic
spatio-temporal process for velocity, and accommodates data assimilation and
uncertainty quantification naturally in the Bayesian inference paradigm.
However, as with many spatio-temporal Bayesian models, this method is fairly
computationally intensive and is somewhat limited in that the latent dynamical
process is assumed to be linear.

We consider here an alternative approach for efficient physics-informed neural
models that provides realistic data-driven estimates of fire front velocities
as well as realistic uncertainty quantification. Specifically, we utilize a
reservoir computing recurrent neural model known as an echo state network
(ESN). Hybrid ESN/statistical approaches for spatio-temporal data were
introduced into statistics by McDermott and Wikle, (2017, 2019) and have
subsequently been used in numerous applications (e.g., Bonas and Castruccio,,
2021; Huang et al.,, 2021; Barredo Arrieta et al.,, 2022). Rather than
estimate neural weight parameters, ESNs and other reservoir models draw
weights randomly from a reservoir distribution, and only train the output
layer parameters. This provides a distinct advantage over traditional neural
models in that much less data is required to train ESNs, making them ideal for
situations such as real-world fire spread where fire boundary data are only
available sporadically (e.g., daily). As with most neural methods, two big
challenges with ESNs relate to hyperparameter tuning (e.g., Wu et al.,, 2018;
Probst et al.,, 2018) and uncertainty quantification (e.g., McDermott and
Wikle,, 2019).

In this paper, we propose an ensemble ESN within the level set framework,
which combines the mechanistically-motivated dynamic model (the level set
method) with reservoir modeling (ESN) of velocities. There are several
contributions of this work. Retaining the level set method as the physically-
motivated dynamical core, our model enables a simple representation of any
complicated boundary and its topological change through time. The ESN, nested
within the level set method to learn propagation speeds, can improve
forecasting performance by accommodating more complicated (i.e., nonlinear)
dynamics that are present in many dynamic systems (Chen and Billings,, 1992).
Importantly, the low computational cost of ESNs makes our model suitable and
attractive for practical use. In that regard, our method addresses the neural
hyperparameter tuning issue efficiently by randomly sampling the
hyperparameters from plausible distributions to construct forecast ensembles.
Finally, our method can quantify the uncertainty in forecasting by using the
ensembles with novel calibration. Compared to forecast intervals derived from
the Bayesian approach of Yoo and Wikle, (2022), our method provides less
formal forecast intervals, but they are effective and much less
computationally expensive. We demonstrate the forecasting performance of the
model on two simulated datasets generated from actual wind fields and two
wildfires in California – the Thomas and Haypress megafires.

The paper is organized as follows. Section 2 provides a brief background on
ESNs and the level set method. Then, in Section 3, our proposed hybrid level
set/ESN model is described, with a novel calibration algorithm to quantify the
uncertainty. We show the forecasting performance of our proposed model through
simulation experiments and real data in Section 4 and 5, respectively.
Finally, Section 6 provides a brief discussion.

##  2 Background

This section provides brief background on ESNs and the level set method to
facilitate discussion of the new model in Section 3.

###  2.1 Echo State Network (ESN)

The ESN has been an appealing type of recurrent neural network due to its
relatively inexpensive computational cost, reduced training data burden, and
conceptual simplicity (Lukoševičious,, 2012). The basic vanilla ESN model for
time $t=1,...,T$ can be written as follows:

$\displaystyle\text{response}:$
$\displaystyle\quad\bm{v_{t}}=g_{o}(\bm{o}_{t})$ $\displaystyle\text{output}:$
$\displaystyle\quad\bm{o}_{t}=\bm{W}^{out}\bm{h}_{t}$
$\displaystyle\text{hidden state}:$
$\displaystyle\quad\bm{h}_{t}=(1-\alpha_{\ell})\bm{h}_{t-1}+\alpha_{\ell}\tilde{\bm{h}}_{t}$
$\displaystyle\quad\tilde{\bm{h}}_{t}=g_{h}\bigg{(}\frac{\nu}{|\lambda_{w}|}\bm{W}^{in}\bm{h}_{t-1}+\bm{U}\bm{x}_{t}\bigg{)},$
(1) $\displaystyle\text{parameters}:$
$\displaystyle\quad\bm{W}^{in}=[w_{i,\ell}^{in}]_{i,\ell}:\gamma_{i,\ell}^{w}Unif(-a_{w},a_{w})+(1-\gamma_{i,\ell}^{w})\delta_{0}$
$\displaystyle\quad\bm{U}=[u_{i,j}]_{i,j}:\gamma_{i,j}^{u}Unif(-a_{u},a_{u})+(1-\gamma_{i,j}^{u})\delta_{0}$
$\displaystyle\quad\gamma_{i,\ell}^{w}\sim Bern(\pi_{w})$
$\displaystyle\quad\gamma_{i,j}^{u}\sim Bern(\pi_{u}),$

where $N\times 1$ vector $\bm{v}_{t}$, $N\times 1$ vector $\bm{o}_{t}$,
$J\times 1$ vector $\bm{h}_{t}$, and $n_{u}\times 1$ vector $\bm{x}_{t}$
correspond to the response, output, hidden units, and input signal,
respectively. Critically, the $J\times J$ matrix $\bm{W}^{in}$ and $J\times
n_{u}$ matrix $\bm{U}$ are pre-specified weight matrices with elements drawn
randomly from specified distributions with hyperparameters $a_{w}$, $a_{u}$,
$\pi_{w}$ and $\pi_{u}$. $\delta_{0}$ is a Dirac delta function at zero. Note
that hyperparameters $\pi_{w}$ and $\pi_{u}$ control the sparsity of the
matrices $\bm{W}^{in}$ and $\bm{U}$, respectively. Additionally, $\lambda_{w}$
is the spectral radius of $\bm{W}^{in}$, $\alpha_{\ell}\in(0,1]$ is the so-
called “leaking rate”, and $\nu$ is a tuning parameter. The element-wise
activation functions $g_{o}(\cdot)$ and $g_{h}(\cdot)$ are user-specified, but
$tanh(\cdot)$ functions are typically chosen for $g_{h}(\cdot)$, an identity
function is typically chosen for $g_{o}(\cdot)$ when the response is
continuous, and a softmax functions is typically chosen for categorical
responses. As discussed by Lukoševičious, (2012), the role of the hidden units
in an ESN is to project the inputs into a higher dimensional space through
nonlinear expansion and to simultaneously remember the input.

To approximate the nonlinear system well, an ESN should exhibit the “echo
state property” (Zhang et al.,, 2012), which essentially ensures that the
reservoir model will eventually not be sensitive to initial conditions. This
is ensured in most cases if the spectral radius of $\bm{W}^{in}$ is less than
1 (Lukoševičious,, 2012). The ratio of the parameter $\nu$ and the largest
eigenvalue $\lambda_{w}$ of $\bm{W}^{in}$ in the hidden state control the
effective spectral radius. That is, $0<\nu<1$ guarantees the spectral radius
of the model will be less than 1.

The key to reservoir methods such as the ESN is that one only needs to
estimate the $N\times J$ read-out matrix $\bm{W}^{out}$, which provides
tremendous savings in terms of computational cost and the amount of data
required to train the model. However, to prevent overfitting, it is very
important that this estimation is regularized, which is typically accomplished
using an $\ell_{2}$ penalty (i.e., ridge regression). For illustration, one
can write the loss function, assuming the identity function for
$g_{o}(\cdot)$, as

$\displaystyle\widehat{\bm{W}}^{out}=\underset{\bm{V}}{\arg\min}\bigg{[}\frac{1}{N}\sum_{t=1}^{T}\bigg{(}(\bm{v_{t}}-\bm{W}^{out}\bm{h_{t}})^{\top}(\bm{v_{t}}-\bm{W}^{out}\bm{h_{t}})\bigg{)}+\tau\sum_{j=1}^{J}||\bm{w}^{out}_{\cdot
j}||^{2}\bigg{]},$ (2)

where $\bm{w}^{out}_{\cdot j}$ is $j^{th}$ column in $\bm{W}^{out}$,
$\tau\in(0,\infty)$ is the regularization parameter, and $||\cdot||$ is the
$\ell_{2}$ norm. One can easily obtain the solution for (2) in closed form
(Lukoševičious,, 2012).

Despite its conceptual simplicity and computational efficiency, selection of
the tuning parameters in the ESN can be burdensome. Lukoševičious, (2012)
provides concise guidelines for the hyperparameters $\pi_{w},\pi_{u},a_{w}$,
and $a_{u}$. In particular, they recommend sparse connections between units
(i.e., small $\pi_{w}$ and $\pi_{u}$) to induce sparsity in $\bm{W}^{in}$ and
$\bm{U}$, but recommend that $\bm{U}$ be more dense than $\bm{W}^{in}$. The
role of $a_{u}$ and $a_{w}$ is to determine the effect of input $\bm{x_{t}}$
and previous hidden state $\bm{h_{t-1}}$ on the current hidden state
$\bm{h_{t}}$ (Lukoševičious,, 2012). In addition to $\pi_{w},\pi_{u},a_{w}$,
and $a_{u}$, one must select the reservoir size $J$, the leakage rate
$\alpha_{\ell}$, and the echo-state parameter $\nu$. Selection of the
hyperparameters $\alpha_{\ell}$ and $J$ can be accomplished by cross-
validation or validation, but can be computationally inefficient. McDermott
and Wikle, (2017) used grid search with a validation approach for finding the
best set of hyperparameters, and McDermott and Wikle, (2019) used a genetic
algorithm (e.g., Sivanandam and Deepa,, 2008) in a “deep” ESN to accommodate
the increased number of hyperparameters present in deep ESNs. Ribeiro et al.,
(2020) applied a Bayesian optimization algorithm to an ESN for short-term load
forecasting, and Thiede and Parlitz, (2019) proposed a gradient-based
optimization algorithm, which was shown to be superior to grid search for
large numbers of hyperparameters. A systematic review of hyperparameter
optimization for ESNs can be found in Matzner, (2022).

Historically, the quantification of uncertainty in neural models presents a
significant challenge (e.g., see Abdar et al.,, 2021; Wikle and Zammit-
Mangion,, 2022, for general discussion). This is also true with reservoir
methods. One useful approach that has been used for ESNs is to generate
ensembles from the reservoir (Yao et al.,, 2013; Sheng et al.,, 2013). For
example, McDermott and Wikle, (2017) and McDermott and Wikle, (2019)
constructed forecasting ensembles by sampling the recurrent and input weight
matrices independently from their reservoir distributions. In a Monte Carlo
sense, the (1-$\alpha$)% forecasting intervals were then obtained from the
element-wise percentiles of the forecasting ensembles. Bonas and Castruccio,
(2021) also generated ensemble forecasts on the training set and quantified
the uncertainty by calibrating the residual with monotone cubic spline
interpolation. In addition to the ensemble approach, McDermott and Wikle,
(2019) adopted stochastic search variable selection priors (George and
McCulloch,, 1993) for their Bayesian ESN. Atencia et al., (2020) employed
Monte Carlo Dropout (Gal and Ghahramani,, 2016) as another straightforward
strategy, which is similar to the ensemble approaches in principle.

###  2.2 The Level Set Method

The level set method (or the level set equation), proposed by Osher and
Sethian, (1988), is an advection-diffusion equation used to model the
evolution of an implicit function in an arbitrary dimension. That is,

$\displaystyle\frac{\partial\phi(\bm{s},t)}{\partial
t}+\bm{v}(\bm{s},t)\cdot\nabla\phi(\bm{s},t)=0,$ (3)

where $\phi(\bm{s},t)$ is an implicit function evaluated at location $\bm{s}$
and time $t$, $\frac{\partial\phi(\bm{s},t)}{\partial t}$ denotes the temporal
partial derivative, and
$\bm{v}(\bm{s},t)=(v_{x}(\bm{s},t),v_{y}(\bm{s},t))^{\top}\in\mathbb{R}^{2}$
is the velocity that controls the direction and speed of the implicit function
evolution. By construction, the level set method requires an implicit
function, leading to the Eulerian representation of the boundary in wildfire
applications. In Eulerian approaches, one represents the boundary implicitly
through a contour function, as opposed to Lagrangian frameworks where one
defines a set of points on the boundary to describe the boundary explicitly.
(Osher and Fedkiw,, 2002). Eulerian representations can be computationally
more expensive than Lagrangian methods, yet Eulerian approaches are typically
more effective in representing complicated boundaries and their change through
time; for example, when wildfires exhibit topological changes such as merging
or splitting.

Figure 1: Illustration of a signed distance function: Three boundaries at
$t=1,2,3$ are depicted from the left to right.

It is common to use a signed distance function for the implicit function in
the level set method due to its straightforward interpretation and
mathematical simplicity. The signed distance function $\phi(\bm{s})$ for
$\forall\bm{s}\in\mathbb{R}^{2}$ is

$\displaystyle\phi(\bm{s})=\begin{cases}-d(\bm{s}),&\quad\bm{s}\in\Omega^{-}\\\
0,&\quad\bm{s}\in\partial\Omega\\\
d(\bm{s}),&\quad\bm{s}\in\partial\Omega^{+}\end{cases},$ (4)

where $\Omega^{-},\Omega^{+}$, and $\partial\Omega$ are the interior, and
exterior of the boundary and the boundary itself, respectively. Here,
$d(\bm{s})$ is the closest distance from $\bm{s}$ to the boundary. Assuming a
distance given by the $\ell_{2}$ norm,

$\displaystyle d(\bm{s})=\min(||\bm{s}-\bm{s}_{I}||),$ (5)

where $\bm{s}_{I}\in\partial\Omega$ is the set of points on the boundary. For
an illustration of a signed distance function evolving in time, see Figure 1.
During the evolution, one can see that a new boundary appears at $t=2$ on the
bottom left, and a closed loop exists inside the main boundary at $t=3$.
Regardless of any topological change the boundary exhibits, one can represent
the boundary with the set $\partial\Omega=\\{\bm{s}:\phi(\bm{s})=0\\}$.
Additionally, a signed distance function has the convenient mathematical
property that

$\displaystyle||\nabla\phi(\bm{s},t)||=1,$ (6)

where $\nabla\phi(\bm{s},t)$ indicates a spatial gradient of $\phi(\bm{s},t)$.

Notably, the frequently-made assumption that the evolution is in the normal
direction to the boundary (e.g., Mallet et al.,, 2009; Alessandri et al.,,
2021; Yoo and Wikle,, 2022) and (6) simplifies (3) so that

$\displaystyle\frac{\partial\phi(\bm{s},t)}{\partial t}+v_{n}(\bm{s},t)$
$\displaystyle=0,$ (7)

where $v_{n}(\bm{s},t)$ is the scalar velocity (speed) in the normal
direction. (7) can then be simply approximated by a forward Euler
discretization

$\displaystyle\frac{\phi(\bm{s},t+\triangle t)-\phi(\bm{s},t)}{\triangle
t}+v_{n}(\bm{s},t)$ $\displaystyle=0,$ (8)

where $\triangle t$ is time-step size. Although one can still solve (3)
without assumption (6) with numerical solvers (Osher and Fedkiw,, 2002), it
induces more computational cost. Importantly, there is not much to be gained
by such complexity in data-driven approaches so long as one can learn the
spatial and temporally varying speed in the normal direction,
$v_{n}(\bm{s},t)$. Indeed, as discussed in Section 1, it is crucial in
wildfire applications to estimate this to obtain realistic fire boundary
propagation in complex fires. Building on the flexibility and convenience of
the ESN, Section 3 presents a novel model that embeds an ESN within the level
set equation to estimate $v_{n}(\bm{s},t)$ in (8).

##  3 Hybrid Level Set/ESN Model

This section presents a model in which an ESN, used to learn propagation
speed, is embedded within the level set method. Section 3.1 illustrates the
structure of the model. Section 3.2, 3.3, and 3.4 present the hyperparameter
tuning approach, uncertainty quantification method, tuning selection, and the
algorithm, respectively.

###  3.1 Nested ESN Model

The hybrid level set/ESN model considers a vectorized form of equation (8)
with additive error in which the spatio-temporal speeds in the normal
direction to the boundary are specified via an ESN. Specifically, for spatial
locations $\bm{s_{i}},1\leq i\leq N$ in the domain of interest
$\mathcal{D}\subset\mathbb{R}^{2}$ and time indexed by $t$, $1\leq t\leq T$,
the full model is given by:

$\displaystyle\text{response}:$
$\displaystyle\quad\bm{\phi}_{t}=\bm{\phi}_{t-\triangle
t}-\bm{v}_{n,t-\triangle t}\triangle t+\bm{\epsilon}_{t}$
$\displaystyle\text{rate of spread}:$ $\displaystyle\quad\bm{v}_{n,t-\triangle
t}=\bm{W}^{out}\bm{h}_{t-\triangle t}$ $\displaystyle\text{hidden state}:$
$\displaystyle\quad\bm{h}_{t-\triangle
t}=(1-\alpha_{\ell})\bm{h}_{t-2\triangle
t}+\alpha_{\ell}\tilde{\bm{h}}_{t-\triangle t}$
$\displaystyle\quad\tilde{\bm{h}}_{t-\triangle
t}=g_{h}\bigg{(}\frac{\nu}{|\lambda_{w}|}\bm{W}^{in}\bm{h}_{t-2\triangle
t}+\bm{U}\bm{x}_{t-\triangle t}\bigg{)}$ (9) $\displaystyle\text{parameters}:$
$\displaystyle\quad\bm{W}^{in}=[w_{i,\ell}^{in}]_{i,\ell}:\gamma_{i,\ell}^{w}Unif(-a_{w},a_{w})+(1-\gamma_{i,\ell}^{w})\delta_{0}$
$\displaystyle\quad\bm{U}=[u_{i,j}]_{i,j}:\gamma_{i,j}^{u}Unif(-a_{u},a_{u})+(1-\gamma_{i,j}^{u})\delta_{0}$
$\displaystyle\quad\gamma_{i,\ell}^{w}\sim Bern(\pi_{w})$
$\displaystyle\quad\gamma_{i,j}^{u}\sim Bern(\pi_{u}),$

where $\bm{\phi}_{t}=(\phi(\bm{s_{1}},t),\ldots,\phi(\bm{s_{N}},t))^{\top}$,
$\bm{v}_{n,t-\triangle t}=(v_{n}(\bm{s_{1}},t-\triangle
t),\ldots,v_{n}(\bm{s_{N}},t-\triangle t))^{\top}$, and
$\bm{\epsilon}_{t}=(\epsilon(\bm{s_{1}},t),\ldots,\epsilon(\bm{s_{N}},t))^{\top}$
denote an observed signed distance function, rate of spread in the normal
direction to the boundary, and error term on the space $\mathcal{D}$ at the
indexed time, respectively. Note that rate of spread (speed), hidden state,
and parameters model are equivalently defined as in (2.1).

The nested ESN model used to estimate the rate of spread
$\bm{v}_{n,t-\triangle t}$ can provide more data-driven learning flexibility
than the linear formulation in Yoo and Wikle, (2022) and the Rothermel model
(Rothermel,, 1972). Importantly, one only needs to estimate the read-out
matrix $\bm{W}^{out}$ after the appropriate hyperparameters are specified. If
ridge regression is used as in (2), one also must specify the penalty
coefficient $\tau$, but this is still computationally efficient.

As indicated in Section 2, uncertainty quantification and the selection of
hyperparameters in this model must be accounted for. We quantify forecast
uncertainty through ensembles, but our method is motivated by the highest
density region proposed by Hyndman, (1996) and calibration strategy of Daw and
Wikle, (2022) as described in Section 3.2. We adopt a random search approach
for hyperparameter tuning (Bergstra and Bengio,, 2012) as described in Section
3.3. The complete algorithm for our hybrid models is presented in Section 3.4.

###  3.2 Uncertainty Quantification

Once the model is fitted, a one-step ahead forecast for time $T+1$ is given by

$\displaystyle\widehat{\bm{\phi}}_{T+1}=\bm{\phi}_{T}-\widehat{\bm{W}}^{out}\bm{h}_{T}\triangle
t,$ (10)

where $\bm{h}_{T}$ is obtained from (3.1) given $\bm{W}^{in}$, $\bm{U}$,
$\alpha_{\ell}$, $J$ and $\nu$. To quantify the uncertainty associated with
this forecast, we take the ensemble approach motivated by Sheng et al.,
(2013),Yao et al., (2013), and McDermott and Wikle, (2017). That is, we
collect an ensemble of $N_{ensemble}$ forecasts
$\widehat{\bm{\phi}}_{en,T+1}=(\widehat{\bm{\phi}}_{T+1}^{(1)},\ldots,\widehat{\bm{\phi}}_{T+1}^{(N_{ensemble})})$
where each forecast is associated with different hyperparameters and reservoir
matrices. We then construct $(1-\alpha)$% forecast intervals at location
$\bm{s}$ as

$\displaystyle\bigg{(}\hat{\ell}_{T+1}(\bm{s}),\hat{u}_{T+1}(\bm{s})\bigg{)}=\underset{\ell_{T+1}(\bm{s}),u_{T+1}(\bm{s})}{\arg\min}\bigg{(}u_{T+1}(\bm{s})-\ell_{T+1}(\bm{s})\bigg{)}$
(11)

such that

$\displaystyle\frac{1}{N_{ensemble}}\sum_{i=1}^{N_{ensemble}}I\bigg{(}\ell_{T+1}(\bm{s})\leq\widehat{\phi}^{(i)}_{T+1}(\bm{s})\leq
u_{T+1}(\bm{s})\bigg{)}=1-\alpha,$ (12)

where $\ell_{T+1}(\bm{s})$ and $u_{T+1}(\bm{s})$ are the interval lower and
upper bounds.

The motivation for (12) is from the calibration strategy of Daw and Wikle,
(2022). To quantify the uncertainty for spatial prediction in reservoir
models, Daw and Wikle, (2022) first obtained the median and inter-quartile
range from ensembles and trained the optimal calibration cutoff value to
ensure a $(1-\alpha)$% coverage rate over the training set. Then, the
prediction intervals for their test set were constructed from the median and
interquartile range from the ensembles over the test set and the cutoff value
obtained from the training set. Here, we do not require a separate training
set calibration because the coverage rate is given directly in (12) based on
ensembles of forecasts.

The highest density region proposed by Hyndman, (1996) motivates the
relationship in (11). The highest density region provides the shortest
intervals, which may consist of several disjoint intervals for multimodal
distributions (Hyndman,, 1996). However, several disjoint intervals complicate
the interpretation in wildfire applications, which leads to quantifying the
uncertainty by single intervals as given in (11). The aim of using the
(1-$\alpha$)% highest density intervals is to exclude the poor forecasts
(e.g., outliers), which we find effective in our simulation experiments and
real data analysis.

###  3.3 Hyperparameter tuning

Hyperparameter selection is one of the most important components when building
deep neural models (Probst et al.,, 2018). The most intuitive and
straightforward approach is to use a grid search to explore the parameter
space. However, such methods are inefficient for relatively high-dimensional
parameter spaces. An efficient alternative is to use random search (e.g.,
Bergstra and Bengio,, 2012), which improves on grid search by randomly
searching over hyperparameters from specified “prior” distributions (Yu and
Zhu,, 2020), and can be implemented at a fraction of the cost of optimization
algorithms. For an extensive review of hyperparameter optimization strategies,
refer to Yu and Zhu, (2020).

We use a version of random search for our hybrid level set/ESN algorithm. The
ensemble approach to ESN modeling allows one to not only sample from the
reservoirs for each ensemble member, but also the hyperparameter
distributions. Specifically, we draw $a_{u}$, $a_{w}$, $J$, $\tau$, $\nu$, and
$\alpha_{\ell}$ for each ensemble member independently from their prior
distributions. The guidelines from Lukoševičious, (2012) can help determine
fixed values for $\pi_{w}$, $\pi_{u}$, and the distributions from which
hyperparameters are drawn. It is important to note that significant
computational efficiency is gained from selecting the hyperparameters through
random search, and also provides uncertainty quantification as described in
Section 3.2.

###  3.4 Hybrid level set/ESN algorithm

Algorithm 1 Hybrid Level Set/ESN Algorithm

1:Input Data $\\{\bm{\phi_{t}}\\}$ for $t=1,...,T$. The prior distributions of
$a_{u}$, $a_{w}$, $J$, $\tau$, $\nu$, and $\alpha_{\ell}$ denoted by
$P(a_{u})$, $P(a_{w})$, $P(J)$, $P(\tau)$, $P(\nu)$, $P(\alpha_{\ell})$,
respectively. Fixed constants $\pi_{w}$ and $\pi_{u}$.

2:for $k=1:N_{ensemble}$ do

3: Randomly draw $a_{u}^{(k)}$, $a_{w}^{(k)}$, $J^{(k)}$, $\tau^{(k)}$,
$\nu^{(k)}$, and $\alpha_{\ell}^{(k)}$ from prior distributions.

4: Construct $\bm{W}^{in,(k)}$, $\bm{U}^{(k)}$ with sampled $a_{w}^{(k)}$,
$a_{u}^{(k)}$, $\nu^{(k)}$, $\pi_{w}$, and $\pi_{u}$.

5: Fit the model in (3.1) with $J^{(k)}$ and $\alpha_{\ell}^{(k)}$, and
estimate $\widehat{\bm{W}}^{out,(k)}$ by ridge regression with
$\lambda^{(k)}$.

6: Obtain the forecasts $\widehat{\bm{\phi}}_{T+1}^{(k)}$ using
$\widehat{\bm{W}}^{out,(k)}$.

7:end for

8:Define the ensemble of forecasts:
$\widehat{\bm{\phi}}_{en,T+1}=\\{\widehat{\bm{\phi}}_{T+1}^{(1)},\ldots,\widehat{\bm{\phi}}_{T+1}^{(N_{ensemble})}\\}$

9:Ensemble forecast median:
$\widehat{\bm{\phi}}_{med,T+1}$=median($\widehat{\bm{\phi}}_{en,T+1}$)

10:Ensemble uncertainty: ($\hat{\ell}_{T+1}$, $\hat{u}_{T+1}$) obtained by
(11) and (12) corresponding to the significance level $\alpha$.

11:output $\widehat{\bm{\phi}}_{med,t+1}$ and ($\hat{\ell}_{T+1}$,
$\hat{u}_{T+1}$)

The hybrid level set/ESN model, including the uncertainty quantification and
hyperparameter selection described in Section 3.2 and 3.3, is detailed in
Algorithm 1.

Step 3 of Algorithm 1 draws hyperparameters $a_{u}^{(k)}$, $a_{w}^{(k)}$,
$J^{(k)}$, $\tau^{(k)}$, $\nu^{(k)}$, and $\alpha_{\ell}^{(k)}$ from their
prior distributions (i.e., random search), from which recurrent and input
weight matrices $\bm{W}^{in,(k)}$ and $\bm{U}^{(k)}$ are constructed in Step
4. Fitting model (3.1) and estimating $\widehat{\bm{W}}^{out,(k)}$ with
$\tau^{(k)}$ is conducted in Step 5, and Step 6 produces forecast ensemble
members

$\displaystyle\widehat{\bm{\phi}}_{T+1}^{(k)}=\bm{\phi}_{T}-\widehat{\bm{W}}^{out,(k)}\bm{h}_{T}^{(k)}\triangle
t,$ (13)

where $\bm{h}_{T}^{(k)}$ is constructed by (3.1) with $\bm{W}^{in,(k)}$,
$\bm{U}^{(k)}$, $\alpha_{\ell}^{(k)}$, $J^{(k)}$ and $\nu^{(k)}$. After
$N_{ensemble}$ iterations of Steps 3 to 6, which creates the ensembles
$\widehat{\bm{\phi}}_{en,T+1}=(\widehat{\bm{\phi}}_{T+1}^{(1)},\ldots,\widehat{\bm{\phi}}_{T+1}^{(N_{ensemble})})$,
Algorithm 1 provides the median of $\widehat{\bm{\phi}}_{en,T+1}$ for the
point forecast estimate and quantifies the uncertainty with $(1-\alpha)$%
intervals of $\widehat{\bm{\phi}}_{en,T+1}$ using (11) and (12).

We show the forecasting performance of this algorithm through simulation
experiments and with wildfire data in Sections 4 and 5, respectively.

Figure 2: Simulation experiments: the evolution of boundaries for each
simulation. In the left figure, the evolution to the northeast is prominent.
In the right figure, evolution is equivalent to the first simulation, but
constraint (solid black line) prevents the spread over it. Figure 3:
Simulation experiments: 95% forecast intervals and forecast median. Figures in
the first and second rows correspond to the first and second simulations. The
solid blue line denotes the true boundary, and the solid red line indicates
the predicted boundary.

##  4 Simulation experiments

To make simulations more realistic, we generate two datasets by evolving a
randomly generated initial closed boundary with the actual wind speeds
obtained from the R package _rWind_ (Fernández-López,, 2022). Specifically, we
first obtain twenty series of wind speeds from 2018-06-25 00:00:00 to
2018-07-04 12:00:00 at every twelve hours over longitude and latitude
[90,150]$\times$[5,40]. Then, the initial closed boundary is evolved from
these wind speeds using equation 8. The first generated dataset (see the left
panel in Figure 2) exhibits an evolution that includes disjoint boundaries
merging into one. The second dataset (see the right panel of Figure 2) is
generated with topological constraints, similar to a spatial area that might
prevent wildfire spread (e.g., rivers, lakes, or oceans). Topological
constraints are depicted in a solid black line for the second dataset. On both
figures, the dashed line indicates the boundary we seek to forecast.

We first obtain a signed distance function on a 31 $\times$ 31 grid on the
domain $[90,150]\times[5,40]$. Then, we leave out data at $t=20$ to evaluate
forecasting performance. That is, we fit our proposed model (3.1) on data from
$t=1,\ldots,19$ and obtain a one-step ahead forecast at $t=20$ by Algorithm 1
with $\triangle t=0.1$. For the hyperparameter distributions, we set
$P(a_{u})=P(a_{w})=(0.1,0.5,1)$,
$P(\alpha_{\ell})=(0.01,0.0621,\ldots,0.9478,1)$,
$P(\nu)=(0.1,0.2,\ldots,0.8,0.9)$, $P(\tau)=(0.001,1.112,\ldots,8.889,10)$,
and $P(J)=(50,$ $\ldots,100)$. The values of $P(J)$ are relatively small to
prevent overfitting, so each ensemble member functions as a relatively weak
learner (McDermott and Wikle,, 2017). Finally, we fix $\pi_{w}=0.3$ and
$\pi_{u}=0.5$ based on the guidelines in Lukoševičious, (2012).

The choice of input is application-specific. McDermott and Wikle, (2017) used
embeddings (lagged values of the data) for the input $\bm{x}_{t}$, which we
follow. Specifically, we set

$\displaystyle\bm{x}_{t}=[\bm{\phi}_{t}^{\top},\bm{\phi}_{t-\triangle
t}^{\top}]^{\top}.$ (14)

Embeddings are motivated by a nonlinear Hamilton-Jacobi equation in the level
set method. Namely, the level set method in (7) becomes a nonlinear Hamilton-
Jacobi equation if $v_{n}(\bm{s},t)$ depends on $\phi(\bm{s},t)$ (Osher and
Sethian,, 1988), which occurs in the ESN if we consider embeddings as in (14).
Finally, we set the initial value of the hidden state
$\bm{h}_{1}=\bm{U}\times\bm{x}_{1+\triangle
t}=\bm{U}\times[\bm{\phi}_{1+\triangle
t}^{\top},\bm{\phi}_{1}^{\top}]^{\top}$.

Figure 3 shows the 95% forecast intervals and median from 3,000 ensembles. One
can see that the 95% forecast intervals are quite narrow, but they mainly
capture the true boundary. The northeast part of the boundary shows the most
variability, which aligns with the evolution pattern depicted in Figure 2.
Moreover, the 95% forecast intervals for the second simulation account for the
topological constraints (although the predicted boundary spreads very slightly
over the topological constraints in the lower bound).

We adopt the threat score (TS) as a summary measure (e.g.,Wikle et al.,
(2019)) to assess the quality of the forecast median. Specifically, the TS is
defined as

$\displaystyle TS=\frac{A_{11}}{A_{11}+A_{10}+A_{01}},$ (15)

where $A_{11}$ is the common area where the model predicts an event that
eventually occurs, $A_{10}$ is the area where the model expects an event to
occur but not eventually occur, and $A_{01}$ is the area where an event
occurs, but the model fails to predict its occurrence. The range of TS summary
measure is from 0 to 1, where a value of 1 indicates a perfect forecast. For a
continuous variable such as a signed distance function, one needs a threshold
parameter $\tau$ to determine the occurrence of an event. Therefore, we set
$\tau=0$ since points on the boundary have a signed distance of 0. That is, if
$\phi(\bm{s},t)\leq 0$, we consider that an event occurred on $\bm{s}$ at $t$.
The TS summary measures for the forecast median of the first and second
simulations are 0.969 and 0.963, demonstrating that our model provides good
forecasting performance.

##  5 Wildfire applications

Figure 4: Fire boundaries for various times during the Haypress fire. Although
the rate of spread was relatively slow until 2017-18-28 02:40:00, it
subsequently accelerated, especially to the west at 2017-09-02 05:45:00.
Figure 5: Forecast medians and 95% intervals for the Haypress fire. The solid
blue lines denote the true boundary, while the solid red lines indicates the
predicted boundary. Most of the variability in forecast intervals is observed
on the left part of the boundary.

In this section we investigate the forecasting performance of our model on two
real datasets - the Haypress and Thomas megafires that occurred in California
in 2017.

###  5.1 The Haypress fire

The Haypress fire, ignited in July 2017 and was contained in January 2018. It
was the largest wildfire among the Orleans Complex fires, which burned over
27,000 acres in Siskiyou County, California. We consider the first 11
observations of the Haypress fire boundary from GeoMac database (Geospatial
Multi-Agency Coordination Group,, 2019) for this analysis. Specifically, we
fit the model (3.1) on the first ten observations, leaving out the eleventh
observation to evaluate forecasting performance. Figure 4 shows the wildfire
propagation from the second observation to the eleventh observation, from
which one can see the rapid rate of spread to the north at 2017-08-30 06:30:00
and to the west at 2017-09-02 05:45:00. The forecasting task is more
challenging than the simulation experiments in that a predominant expansion to
the west shown at 2017-09-02 05:45:00 is not observed from 2017-08-13 23:32:00
to 2017-09-01 06:30:00. At 2017-09-01 06:30:00, the Haypress fire shows a
slight spread to the west, but it is not as striking as observed in the
boundary at 2017-09-02 05:45:00.

With a signed distance function on a $31\times 31$ grid over
[-123.5053,-123.3252] $\times$ [41.3840, 41.5471] by (4), we use Algorithm 1
to fit the model (3.1) with the same hyperparameter distributions and
$\triangle t=0.1$ as in the simulation experiments. Also, the inputs and
initial value of the hidden state are chosen in the same manner. Figure 5
shows the 95% forecast intervals and median from 3,000 ensembles. The forecast
median captures the evolution to the west, and the 95% intervals provide good
coverage of the true boundary. Given the small number of observations used to
train the model, the model visually shows very good forecasting performance.

To more objectively compare the forecast performance of our model to that by
Yoo and Wikle, (2022), we measure the quality of forecast intervals using the
interval score (IS) (Gneiting and Raftery,, 2012), which is given by

$\displaystyle S_{\alpha}^{Int}$
$\displaystyle(\ell(\bm{s},t),u(\bm{s},t),\phi(\bm{s},t))=$
$\displaystyle(u(\bm{s},t)-\ell(\bm{s},t))+\frac{2}{\alpha}(\ell(\bm{s},t)-\phi(\bm{s},t))I(\phi(\bm{s},t)<\ell(\bm{s},t))$
(16)
$\displaystyle+\frac{2}{\alpha}(\phi(\bm{s},t)-u(\bm{s},t))I(\phi(\bm{s},t)>u(\bm{s},t)),$

where $\ell(\bm{s},t)$, $u(\bm{s},t)$, $\alpha$, and $\phi(\bm{s},t)$ are the
lower bound, upper bound, significance level, and true signed distance
function on $\bm{s}$ at $t$. The IS consists of the length of an interval
$(u(\bm{s},t)-\ell(\bm{s},t))$ and the penalty
$(\frac{2}{\alpha}(\ell(\bm{s},t)-\phi(\bm{s},t))$ or
$\frac{2}{\alpha}(\phi(\bm{s},t)-u(\bm{s},t)))$ induced if forecast intervals
fail to contain $\phi(\bm{s},t)$. A smaller IS indicates better forecast
intervals. The IS is intuitively appealing since it favors narrow forecast
intervals but penalizes if those are too narrow. Furthermore, the IS is a
proper scoring rule (Gneiting and Raftery,, 2012). We summarize the quality of
forecasting intervals by

$\displaystyle IS=\frac{\sum_{\bm{s}:\phi(\bm{s},t)\leq
0}S_{\alpha}^{Int}(\ell(\bm{s},t),u(\bm{s},t),\phi(\bm{s},t))}{|\\{\bm{s}:\phi(\bm{s},t)\leq
0\\}|},$ (17)

where $|\cdot|$ denotes the cardinality of the set and $t$=2017-09-02 05:45:00
for the Haypress fire. We only consider the locations inside the true boundary
at 2017-09-02 05:45:00 since our primary interest is the boundary (indirectly,
the burned and burning area), not the whole domain.

As given in Table 1, the TS summary measure and IS from model applied to the
Haypress fire forcast are 0.849 and 0.776, respectively, while those from Yoo
and Wikle, (2022) are 0.791 and 1.730, demonstrating better forecast
performance of the hybrid level set/ESN model. We also note that the
computational cost of the hybrid level set/ESN is much lower than Yoo and
Wikle, (2022). Specifically, the algorithm presented in this paper takes less
than 3 minutes on a standard laptop computer compared to 3 hours for the model
in Yoo and Wikle, (2022) on a high-performance cluster.

Table 1: Comparison of Interval Score (IS) and Threat Score (TS) between models. $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ refer to the model presented here and the model by Yoo and Wikle, (2022), respectively. Model |  | The Haypress fire |  | The Thomas fire
---|---|---|---|---
| Interval Score | Threat Score |  | Interval Score | Threat Score
$\mathcal{M}_{1}$ |  | 0.776 | 0.849 |  | 4.901 | 0.964
$\mathcal{M}_{2}$ |  | 1.730 | 0.791 |  | 5.286 | 0.954

###  5.2 The Thomas fire

Figure 6: Thomas fire boundaries for 20 times Thomas fire (note, the first
observation at 2017-12-05 04:40:00 is omitted). Compared to the Haypress fire,
the rate of spread is relatively consistent in general.

The massive Thomas fire in Southern California burned 114,000 hectares from
the time of its ignition on December 4, 2017 to its containment in January
2018. Even after containment, post-fire debris flow continued to present
hazards (Addison and Oommen,, 2020). We consider 21 observations of fire
boundaries from the GeoMac database (Geospatial Multi-Agency Coordination
Group,, 2019) as depicted in Figure 6. Notably, the spread rate is relatively
more uniform than the Haypress fire, especially after 2017-12-10 23:21:00.
However, there was a sudden extension of the boundary at 2017-12-24 18:00:00.
This makes forecasting at 2017-12-24 18:00:00 challenging as the boundaries
between 2017-12-19 00:23:00 and 2017-12-23 23:10:00 were evolving quite
slowly.

We use the first twenty observations to fit the model and evaluate forecast
performance at 2017-12:24 18:00:00. A signed distance function is constructed
on a $30\times 30$ grid over [-119.8,-118.8] $\times$ [34.2, 34.7] by (4).
Using the same hyperparameter distributions, inputs, $\triangle t=0.1$, and
initial value of the hidden state as the simulation experiments, we obtain the
forecast median and 95% forecast intervals from 3,000 ensembles, as shown in
Figure 7. The forecast median fails to predict the upper part of the boundary
connected at 2017-12-24 18:00:00, but importantly, the 95% forecast intervals
capture this feature with a loop in the lower bound. This is a significant
benefit of uncertainty quantification approach presented here.

Figure 7: Forecast median and 95% intervals for the Thomas fire. The solid
blue line denotes the true boundary and the solid red line indicates predicted
boundary

We also compare objectively the hybrid level set/ESN model to the approach in
Yoo and Wikle, (2022) with the TS summary measure and IS. As given in Table 1,
the hybrid model gives TS and IS values of 0.964 and 4.901, respectively,
while Yoo and Wikle, (2022) gave 0.954 and 5.286. These summary measures
demonstrate superiority of the hybrid model, but the difference is not as
dramatic as with the Haypress fire. In addition, we observe visually that the
95% credible intervals from Yoo and Wikle, (2022) cover the true boundary
without a closed loop in the lower bound, which is found in Figure 7.
Nevertheless, the wider credible intervals by Yoo and Wikle, (2022) have more
negative impacts on the IS than the penalty induced by not fully covering the
true boundary in the hybrid model.

##  6 Discussion

Although the level set method with a signed distance function has become a
powerful tool in wildfire applications, the importance of data-driven learning
is often neglected. Given the inhomogeneous environmental conditions and rapid
spread of wildfires, it is important that the forecasting model adapt through
data-driven learning. In addition, an important aspect of a trustworthy
forecasting model is realistic uncertainty quantification, which is not
available with deterministic model applications. Finally, a real-world
wildfire forecasting model should account for nonlinearity in the fire
evolution to accommodate more realistic dynamics.

We formulate our model by combining the mechanistically-motivated dynamic
model (the level set method) with reservoir modeling of fire front propagation
speed (ESN). Importantly, this model allows for data-driven learning as well
as accounting for nonlinearity. Also, our algorithm is highly efficient and
provides realistic forecasting intervals without requiring significant
training costs and large amounts of data. Through simulation experiments, we
demonstrate outstanding forecasting performance. We also show that our model
can provide better forecasting intervals than a state-of-the-art Bayesian
level set model (Yoo and Wikle,, 2022) on the Haypress and Thomas fire.

Note that we used the same distributions for hyperparameters for the
simulations and real-world examples. Although for certain applications one
might need to choose different distributions depending on the scale of data,
our proposed random search algorithm is computationally efficient and relieves
much of the burden in the tuning parameter selection as well as providing
realistic forecasting intervals. Also, note that the same time-step size
$\triangle t$ is used for both numerical experiments and wildfire
applications. We observed that the results are not very sensitive to the
choice of $\triangle t$ as the estimated readout matrix
$\widehat{\bm{W}}^{out}$ adjusts to different values of $\triangle t$. Note,
however, that our hybrid level set model is not limited to wildfire
applications. The methodology presented here can also be used to infer
shrinking boundaries, such as the response to a tumor after treatment, or the
retreat of Antarctic glaciers. The key is the data-driven estimation of the
speed in the normal direction to the boundary, which can easily accommodate
expansion or contraction (or both in different parts of the domain).

Our model shows superior results and is more computational efficient compared
to Yoo and Wikle, (2022), but this does not necessarily mean that our model is
more advantageous than Yoo and Wikle, (2022) in every respect. Specifically,
one limitation of our model is the failure to explicitly determine the effects
of covariates (e.g., slope and landcover). Note that covariates, such as
slope, aspect, existing vegetation cover, and forest canopy cover, can be
incorporated in $\bm{h}_{1}$ as:

$\displaystyle\bm{h}_{1}=\bm{U}\bm{x}_{1+\triangle
t}+\bm{U}_{2}\bm{z_{slp}}+\bm{U}_{3}\bm{z_{asp}}+\bm{U}_{4}\bm{z_{veg}}+\bm{U}_{5}\bm{z_{cnp}},$
(18)

where $\bm{Z}=[\bm{z}_{slp},\bm{z}_{asp},\bm{z}_{veg},\bm{z}_{cnp}]$ denote
covariates for slope, aspect, existing vegetation cover, and forecast canopy
cover, while $\bm{U}_{2},\bm{U}_{3},\bm{U}_{4},\bm{U}_{5}$ are pre-specified
(random) weight matrices. When considering such a model, we did not observe
significant improvement in forecasting with (18) and one cannot perform
inference on the importance of these effects directly, as is possible in the
fully hierarchical Bayesian framework in Yoo and Wikle, (2022). Alternatively,
one can add an additional stage in (3.1) to concatenate covariates in hidden
states as:

$\displaystyle\bm{h}_{t-\triangle t}^{*}=[\bm{h}_{t-\triangle
t}^{\top},\bm{z}^{\top}]^{\top},$ (19)

where
$\bm{z}^{\top}=[z_{slp}^{\top},z_{asp}^{\top},z_{veg}^{\top},z_{cnp}^{\top}]^{\top}$.
Even though the readout matrix corresponding to $(19)$ can quantify the
element-wise effect of covariates, this model did not show significant
improvement in forecasting. Therefore, a subject of future research is to
consider alternative approaches for introducing covariates into our hybrid
model in a manner that can provide explainability regarding covariate
importance (e.g., Wikle et al.,, 2022).

Another direction in which to extend our model is to account for topological
constraints explicitly. Although we demonstrate in Figure 3 that our model
effectively accounts for such constraints implicitly when there is data to
inform the normal speed, in situations where there are data gaps in time, or
multiple time step forecasts are of interest, then explicit constraints would
be important to produce realistic forecasts (e.g., accounting for water
bodies). Employing a re-initialization technique in the model could be another
path. Osher and Sethian, (1988) suggested occasional re-initialization to
maintain a signed distance function. Nevertheless, we observe that re-
initialization is not necessary for one-step ahead forecasting given data at
each time step (see Figure 3, 5, and 7). Yet, for multiple-step ahead
forecasting, re-initialization would likely be beneficial.

Lastly, it would be instructive to combine the level set method with other
variants of ESNs. Our model is based on the basic ESN, but we could also
consider quadratic or deep ESNs (e.g., McDermott and Wikle,, 2017, 2019).

## References

  * Abdar et al., (2021) Abdar, M., Pourpanah, F., Hussain, S., Rezazadegan, D., Liu, L., Ghavamzadeh, M., Fieguth, P., Cao, X., Khosravi, A., Acharya, U., Makarenkov, V., and Nahavandi, S. (2021).  A review of uncertainty quantification in deep learning: Techniques, applications and challenges.  Information Fusion, pages 243–297.  https://doi.org/10.1016/j.inffus.2021.05.008.
  * Addison and Oommen, (2020) Addison, P. and Oommen, T. (2020).  Post-fire debris flow modeling analyses: case study of the post-thomas fire event in california.  Natural Hazards, 100:329–343.  https://doi.org/10.1007/s11069-019-03814-x.
  * Alessandri et al., (2021) Alessandri, A., Bagnerini, P., Gaggero, M., and Mantelli, L. (2021).  Parameter estimation of fire propagation models using level set methods.  Applied Mathematical Modelling, 92:731–747.  https://doi.org/10.1016/j.apm.2020.11.030.
  * Atencia et al., (2020) Atencia, M., Stoean, R., and Joya, G. (2020).  Uncertainty quantification through dropout in time series prediction by echo state networks.  Mathematics, 8:1374.  https://doi.org/10.3390/math8081374.
  * Banerjee et al., (2015) Banerjee, S., Carlin, B., and Gelfand, A. (2015).  Hierarchical Modeling and Analysis for Spatial Data.  Chapman & Hall/CRC, Boca Raton, FL, USA, 2 edition.  https://doi.org/10.1201/b17115.
  * Barredo Arrieta et al., (2022) Barredo Arrieta, A., Gil-Lopez, S., Laña, I., Bilbao, M. N., and Del Ser, J. (2022).  On the post-hoc explainability of deep echo state networks for time series forecasting, image and video classification.  Neural Computing and Applications, 34(13):10257–10277.  https://doi.org/10.1007/s00521-021-06359-y.
  * Bergstra and Bengio, (2012) Bergstra, J. and Bengio, Y. (2012).  Random search for hyper-parameter optimization.  Journal of Machine Learning Research, 13:281–305.
  * Bonas and Castruccio, (2021) Bonas, M. and Castruccio, S. (2021).  Calibration of spatial forecasts from citizen science urban air pollution data with sparse recurrent neural networks.  arXiv preprint arXiv:2105.02971.
  * Brotons et al., (2013) Brotons, L., Aquilué, N., de Cáceres, M., Fortin, M.-J., and Fall, A. (2013).  How fire history, fire suppression practices and climate change affect wildfire regimes in mediterranean landscapes.  PloS one, 8(5):e62392.  https://doi.org/10.1371/journal.pone.0062392.
  * Burge et al., (2022) Burge, J., Bonanni, M. R., Hu, R., and Ihme, M. (2022).  Recurrent convolutional deep neural networks for modeling time-resolved wildfire spread behavior.  arXiv preprint arXiv:2210.16411.
  * Cannon and DeGraff, (2009) Cannon, S. and DeGraff, J. (2009).  The increasing wildfire and post-fire debris-flow threat in western usa, and implications for consequences of climate change.  In Landslides – Disaster Risk Reduction, pages 177–190. Springer.  https://doi.org/10.1007/978-3-540-69970-5_9.
  * Chen and Billings, (1992) Chen, S. and Billings, S. (1992).  Neural networks for nonlinear dynamic system modelling and identification.  International Journal of Control, 56(2):319–346.  https://doi.org/10.1080/00207179208934317.
  * Cressie and Wikle, (2011) Cressie, N. and Wikle, C. K. (2011).  Statistics for Spatio-Temporal Data.  Wiley, New Jersey.
  * (14) Dabrowski, J. J., Huston, C., Hilton, J., Mangeon, S., and Kuhnert, P. (2022a).  Towards data assimilation in level-set wildfire models using bayesian filtering.  arXiv preprint arXiv:2206.08501.
  * (15) Dabrowski, J. J., Pagendam, D. E., Hilton, J., Sanderson, C., MacKinlay, D., Huston, C., Bolt, A., and Kuhnert, P. (2022b).  Bayesian physics informed neural networks for data assimilation and spatio-temporal modelling of wildfires.  arXiv preprint arXiv:2212.00970.
  * Daw and Wikle, (2022) Daw, R. and Wikle, C. K. (2022).  Reds: Random ensemble deep spatial prediction.  Environmetrics, page e2780.  https://doi.org/10.1016/j.spasta.2020.100408.
  * Fendell and Wolff, (2001) Fendell, F. and Wolff, M. (2001).  Wind-aided fire spread.  In Forest Fires: Behavior and Ecological Effects, page 171–223. Academic Press.  10.1016/B978-012386660-8/50008-8.
  * Fernández-López, (2022) Fernández-López, J. (2022).  rwind: Download, edit and include wind and sea currents data in ecological and evolutionary analysis.  R package version 1.1.7.
  * Gal and Ghahramani, (2016) Gal, Y. and Ghahramani, Z. (2016).  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  arXiv preprint arXiv:1506.02142.
  * George and McCulloch, (1993) George, E. and McCulloch, R. (1993).  Variable selection via gibbs sampling.  Journal of the American Statistical Association, 88:881–889.  https://doi.org/10.2307/2290777.
  * Geospatial Multi-Agency Coordination Group, (2019) Geospatial Multi-Agency Coordination Group (2019).  Historic perimeters combined 2000-2018 geomac.  Geospatial Multi-Agency Coordination Group.
  * Gibou et al., (2018) Gibou, F., Fedkiw, R., and Osher, S. (2018).  A review of level-set methods and some recent applications.  Journal of Computational Physics, 353:82–109.  https://doi.org/10.1016/j.jcp.2017.10.006.
  * Gneiting and Raftery, (2012) Gneiting, T. and Raftery, A. (2012).  Strictly proper scoring rules, prediction, and estimation.  Journal of the American Statistical Association, 102:359–378.  https://doi.org/10.1198/016214506000001437.
  * Hohner et al., (2019) Hohner, A., Rhoades, C., Wilkerson, P., and Rosario-Ortiz, F. (2019).  Wildfires alter forest watersheds and threaten drinking water quality.  Accounts of Chemical Research, 52:1234–1244.  DOI:10.1021/acs.accounts.8b00670.
  * Huang et al., (2021) Huang, H., Castruccio, S., and Genton, M. G. (2021).  Forecasting high-frequency spatio-temporal wind power with dimensionally reduced echo state networks.  arXiv preprint arXiv:2102.01141.
  * Hyndman, (1996) Hyndman, R. J. (1996).  Computing and graphing highest density regions.  The American Statistician, 50:120–126.  https://doi.org/10.2307/2684423.
  * Lo, (2012) Lo, S.-E. (2012).  A fire simulation model for heterogeneous environments using the level set method.  Claremont Graduate University, 2012 Ph.D. thesis.  10.5642/cguetd/72.
  * Lukoševičious, (2012) Lukoševičious, M. (2012).  A practical guide to applying echo state networks.  In Neural Networks: Tricks of the Trade, pages 659–686. Springer.  https://doi.org/10.1007/978-3-642-35289-8_36.
  * Mallet et al., (2009) Mallet, V., Keyes, D. E., and Fendell, F. E. (2009).  Modeling wildland fire propagation with level set methods.  Computers and Mathematics with Applications, 57(7):1089–1101.  https://doi.org/10.1016/j.camwa.2008.10.089.
  * Matzner, (2022) Matzner, F. (2022).  Hyperparameter tuning in echo state networks.  arXiv preprint arXiv:2207.07976.
  * McDermott and Wikle, (2017) McDermott, P. and Wikle, C. (2017).  An ensemble quadratic echo state network for non-linear spatio-temporal forecasting.  Stat, 6:315–330.  https://doi.org/10.1002/sta4.160.
  * McDermott and Wikle, (2019) McDermott, P. and Wikle, C. K. (2019).  Deep echo state networks with uncertainty quantification for spatio-temporal forecasting.  Environmetrics, 30:e2553.  https://doi.org/10.1002/env.2553.
  * Miller et al., (2015) Miller, C., Hilton, J., Sullivan, A., and Prakash, M. (2015).  Spark–a bushfire spread prediction tool.  In International Symposium on Environmental Software Systems, pages 262–271. Springer.  10.1007/978-3-319-15994-2_26.
  * Osher and Fedkiw, (2002) Osher, S. and Fedkiw, R. (2002).  Level Set Methods and Dynamic Implicit Surfaces.  Springer-Verlag, New York.  https://doi.org/10.1007/b98879.
  * Osher and Sethian, (1988) Osher, S. and Sethian, J. A. (1988).  Fronts propagating with curvature dependent speed: Algorithms based on hamilton-jacobi formulations.  Journal of Computational Physics, 79(1):12–49.  https://doi.org/10.1016/0021-9991(88)90002-2.
  * Probst et al., (2018) Probst, P., Boulesteix, A.-L., and Bischl, B. (2018).  Tunability: Importance of hyperparameters of machine learning algorithms.  arXiv preprint arXiv:1802.09596.
  * Raissi et al., (2019) Raissi, M., Perdikaris, P., and Karniadakis, G. E. (2019).  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational physics, 378:686–707.  https://doi.org/10.1016/j.jcp.2018.10.045.
  * Ribeiro et al., (2020) Ribeiro, G., Sauer, J., Fraccanabbia, N., Mariani, V., and Coelho, L. (2020).  Bayesian optimized echo state network applied to short-term load forecasting.  Energies, 13:2390.  https://doi.org/10.3390/en13092390.
  * Rothermel, (1972) Rothermel, R. C. (1972).  A mathematical model for predicting fire spread in wildland fuels.  Res. Pap. INT-115. Ogden, UT: U.S. Department of Agriculture, Intermountain Forest and Range Experiment Station.
  * Sheng et al., (2013) Sheng, C., Zhao, J., Wang, W., and Leung, H. (2013).  Prediction intervals for a noisy nonlinear time series based on a bootstrapping reservoir computing network ensemble.  IEEE Transactions on Neural Networks and Learning Systems, 24(7):1036–1048.  10.1109/TNNLS.2013.2250299.
  * Sivanandam and Deepa, (2008) Sivanandam, S. and Deepa, S. (2008).  Introduction to Genetic Algorithms.  Springer, Berlin.  https://doi.org/10.1007/978-3-540-73190-0.
  * Thiede and Parlitz, (2019) Thiede, L. and Parlitz, U. (2019).  Gradient based hyperparameter optimization in echo state networks.  Neural Networks, 115:23–29.  https://doi.org/10.1016/j.neunet.2019.02.001.
  * Wikle and Zammit-Mangion, (2022) Wikle, C. and Zammit-Mangion, A. (2022).  Statistical deep learning for spatial and spatio-temporal data.  arXiv preprint arXiv:2206.02218.
  * Wikle et al., (2022) Wikle, C. K., Datta, A., Hari, B. V., Boone, E. L., Sahoo, I., Kavila, I., Castruccio, S., Simmons, S. J., Burr, W. S., and Chang, W. (2022).  An illustration of model agnostic explainability methods applied to environmental data.  Environmetrics, page e2772.
  * Wikle et al., (2019) Wikle, C. K., Zammit-Mangion, A., and Cressie, N. (2019).  Spatio-Temporal Statistics with R.  Chapman & Hall/CRC, Florida.
  * Wu et al., (2018) Wu, Q., Fokoue, E., and Kudithipudi, D. (2018).  On the statistical challenges of echo state networks and some potential remedies.  arXiv preprint arXiv:1802.07369.
  * Yao et al., (2013) Yao, W., Zeng, Z., Lian, C., and Tang, H. (2013).  Ensembles of echo state networks for time series prediction.  In 2013 Sixth International Conference on Advanced Computational Intelligence (ICACI), pages 299–304.  10.1109/ICACI.2013.6748520.
  * Yoo and Wikle, (2022) Yoo, M. and Wikle, C. K. (2022).  A bayesian spatio-temporal level set dynamic model and application to fire front propagation.  arXiv preprint arXiv:2210.14978.
  * Yu and Zhu, (2020) Yu, T. and Zhu, H. (2020).  Hyper-parameter optimization: A review of algorithms and applications.  arXiv preprint arXiv:2003.05689.
  * Zhang et al., (2012) Zhang, B., Miller, D. J., and Wang, Y. (2012).  Nonlinear system modeling with random matrices: Echo state networks revisited.  IEEE Transactions on Neural Networks and Learning Systems, 23(1):175–182.  10.1109/TNNLS.2011.2178562.
  * Zuo et al., (2015) Zuo, Z., Shuai, B., Wang, G., Liu, X., Wang, X., Wang, B., and Chen, Y. (2015).  Convolutional recurrent neural networks: Learning spatial dependencies for image representation.  In 2015 IEEE Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), pages 18–26.  10.1109/CVPRW.2015.7301268.


# Blast in a One-Dimensional Cold Gas: From Newtonian Dynamics to
Hydrodynamics

Subhadip Chakraborti1, Santhosh Ganapa1, P. L. Krapivsky2,3 and Abhishek Dhar1
1International Centre for Theoretical Sciences, Tata Institute of Fundamental
Research, Bengaluru 560089, India  2Department of Physics, Boston University,
Boston, Massachusetts 02215, USA  3Skolkovo Institute of Science and
Technology, 143026 Moscow, Russia

###### Abstract

A gas composed of a large number of atoms evolving according to Newtonian
dynamics is often described by continuum hydrodynamics. Proving this
rigorously is an outstanding open problem, and precise numerical
demonstrations of the equivalence of the hydrodynamic and microscopic
descriptions are rare. We test this equivalence in the context of the
evolution of a blast wave, a problem that is expected to be at the limit where
hydrodynamics could work. We study a one-dimensional gas at rest with
instantaneous localized release of energy for which the hydrodynamic Euler
equations admit a self-similar scaling solution. Our microscopic model
consists of hard point particles with alternating masses, which is a
nonintegrable system with strong mixing dynamics. Our extensive microscopic
simulations find a remarkable agreement with Euler hydrodynamics, with
deviations in a small core region that are understood as arising due to heat
conduction.

The Navier-Stokes-Fourier (NSF) equations zeytounian2012  describing the
evolution of the density, velocity, and temperature apply to an enormous range
of phenomena, e.g., to atmospheric flows. At the fundamental level, however,
molecules follow Newton’s equations of motion. How accurate is the
hydrodynamic description and what are the limits of its applicability? Here we
address this question, which is of fundamental and practical importance. A
rigorous derivation of the continuum hydrodynamics description from the
atomistic one is an open problem esposito2004 ; gallagher2019 ; gorban2018
though there are phenomenological derivations via coarse-graining procedures
resibois1977 ; espanol2009 ; spohn2012 . Much progress has been made using the
Boltzmann equation – one derives the equations of hydrodynamics through a
systematic expansion in a small parameter chapman1990 ; de1989 ; saint2009 .
This still leaves one with the problem of deriving the Boltzmann equation from
Newton’s equations which has not been fully achieved even for dilute gases.
For the case where one adds a weak noise to the Newtonian dynamics (still
satisfying the same conservation laws) a rigorous derivation of the Euler
equations for the hydrodynamic fields has been achieved olla1993 .

Surprisingly, there appears to be no direct numerical verification that
hydrodynamics accurately reproduces the predictions of the microscopic
dynamics. The present work provides such a detailed comparison in the context
of the classic blast-wave problem for which a self-similar scaling solution of
the Euler equations was obtained more than sixty years back by Taylor
Taylor19501 ; Taylor19502 , von Neumann VonNeumann1963  and Sedov Sedov1946 ;
Sedov2014 , and is referred to as the TvNS solution. The evolution of a blast
wave emanating from an intense explosion was first studied to understand the
mechanical effect of bombs. The rapid release of a large amount of energy in a
localized region produces a surface of discontinuity beyond which the
quantities concerned like density, velocity, and temperature fields change
discontinuously LandauBook ; ZeldovichBook . The blast wave problem thus
presents an extreme case to test the validity of hydrodynamics. From the point
of microscopic models, the hard sphere system would be the natural candidate
since much is known analytically. Second, this system can be simulated very
efficiently using event-driven simulations. However, large scale molecular
dynamics (MD) simulations Antal2008 ; Jabeen2010 ; Barbier2015 ; Barbier2015a
; Barbier2016 ; Joy2017 ; Joy2019 ; Joy2021  have so far not found clear
agreement with the TvNS solution. It was suggested that possible reasons for
the differences could be the lack of local equilibration or due to the
contribution of viscosity and heat conduction not included in the TvNS
analysis.

In this Letter, we address this question by studying hard point particles with
binary mass distribution — in one dimension, particles with equal masses just
exchange velocities, so mass dispersion is necessary for relaxation, and the
binary mass distribution is the simplest setting where relaxation is possible.
Furthermore, we assume that adjacent particles have different masses (say
$m_{1}$ and $m_{2}$). This alternating hard particle (AHP) gas has been
extensively investigated in the context of the breakdown of Fourier’s law of
heat conduction in one dimension Garrido2001 ; Dhar2001 ; Grassberger2002 ;
Casati2003 ; Cipriani2005 ; Chen2014 ; Hurtado2016 ; Zhao2018 ; Lepri2020 .
The hard particle (and hard rod) system was investigated earlier in the
context of the breakdown of the hydrodynamic description in one dimension
Kadanoff1995 ; Hurtado2006  and more recently the evolution of the AHP
starting from a domain wall initial condition was studied in Mendl2017  and
incomplete thermalization of hard rods in a harmonic trap was observed in
moore2018 . Compared to the hard sphere system in higher dimensions, the 1D
gas of point particles has several advantages – the equation of state is
exactly that of an ideal gas, and simulations are faster since collisions
occur only with nearest neighbors. Note that while the equilibrium physics is
that of an ideal gas, the dynamics is nonintegrable and known to have good
ergodic properties casati1999 .

Figure 1: Heat maps showing the spatiotemporal evolution of the density,
velocity and temperature fields, starting from initial conditions
corresponding to a Gaussian initial temperature profile and
$\rho(x,0)=\rho_{\infty}=1.5,v(x,0)=0$. The simulation parameters were
$N=L=24000$, $E=32,\mu=1.5$, and an ensemble average was taken over $10^{4}$
initial conditions.

We study the evolution of the blast wave initial condition in the AHP gas.
From extensive molecular dynamics simulations of the AHP gas, we compute the
evolution of the density, velocity, and temperature fields and thereby extract
the scaling forms obtained in the long time limit. We make comparisons with
the TvNS scaling solution which we obtain exactly. We find that a complete
explanation of the simulation of the blast requires us to go beyond the Euler
equation and include the effect of heat conduction. We thus discuss the NSF
equations, for which we present results from a numerical solution as well as a
scaling analysis.

The TvNS solution. – We take a 1D ideal gas at zero temperature and uniform
density $\rho_{\infty}$, and suddenly inject energy $E$ into a localized
region of extent $\sigma$. The Euler equations for the density $\rho(x,t)$,
velocity $v(x,t)$, and temperature $T(x,t)$ read

$\displaystyle\partial_{t}\rho+\partial_{x}(\rho v)=0,$ (1a)
$\displaystyle\partial_{t}(\rho v)+\partial_{x}(\rho v^{2}+P)=0,$ (1b)
$\displaystyle\partial_{t}(\rho e)+\partial_{x}(\rho ve+Pv)=0.$ (1c)

For an ideal gas $e=v^{2}/2+T/(2\mu)$ and $P=\rho k_{B}T/\mu$, with
$\mu=(m_{1}+m_{2})/2$ for our binary mass system. We set Boltzmann’s constant
to unity: $k_{B}=1$. A shock wave is formed, and it advances as
$R\sim(Et^{2}/\rho_{\infty})^{1/(d+2)}$ in $d$ dimensions LandauBook , so in
our 1D setting

$R(t)=\left(\frac{Et^{2}}{A\rho_{\infty}}\right)^{\frac{1}{3}},$ (2)

where $A$ is a dimensionless constant factor which, quite remarkably, can be
computed exactly supplemental . From dimensional analysis one can further show
LandauBook  that the fields will take the following self-similar scaling form

$\displaystyle\rho(x,t)$ $\displaystyle=\rho_{\infty}G(\xi),$ (3a)
$\displaystyle v(x,t)$
$\displaystyle=\frac{2}{3}\frac{x}{t}V(\xi)=\frac{2\alpha}{3t^{1/3}}\xi
V(\xi),$ (3b) $\displaystyle T(x,t)$
$\displaystyle=\frac{4\mu}{27}\frac{x^{2}}{t^{2}}Z(\xi)=\frac{4\mu\alpha^{2}}{27t^{2/3}}\xi^{2}Z(\xi).$
(3c)

Here $\xi=x/R(t)$ is the scaled spatial coordinate,
$\alpha=[E/(A\rho_{\infty})]^{1/3}$ and the scaling functions $G,V,Z$ need to
be determined. The factors $2/3,4/27$ are inserted for convenience; e.g., from
Eq. (2) one finds that the velocity of the shock wave is $({2}/{3}){R}/{t}$
and this suggests the use of the factor $2/3$. Plugging these into Eqs.
(1a–1c) we find that the scaling functions satisfy a set of coupled first
order ordinary differential equations in the variable $\xi$. Using the
condition of conservation of energy and the so-called Rankine-Hugoniot
conditions, specifying the field discontinuities at the shock front, allows
one to obtain a complete closed-form solution of the problem, i.e., the
functions $G,V,Z$ and the constant $A$ BarenblattBook . In the Supplemental
Material supplemental  and Ref. arxiv  we describe some details of the
solution.

Microscopic dynamics and initial conditions. –  Our system is a 1D gas of $N$
hard point particles moving inside a box $(-L/2,L/2)$. The only interactions
between particles are through collisions between nearest neighbors that
conserve energy and momentum and also the ordering of the particles. Between
collisions, the particles move ballistically with constant speeds while the
postcollision velocities follow from Newton’s laws (see Ref. supplemental ).
For the AHP gas, the only conserved quantities are particle number, total
momentum, and energy and we expect a hydrodynamic description in terms of the
corresponding conserved fields that are obtained from the microscopic
variables using the standard relations:

$\displaystyle(\rho,p,E)=\sum_{j=1}^{N}\left\langle
m_{j}(1,v_{j},v_{j}^{2}/2)~{}\delta[q_{j}(t)-x]\right\rangle,$ (4)

where $\langle...\rangle$ indicates an average over an initial distribution of
microstates that correspond to the same initial macrostate. We define
$v(x,t)=p/\rho$ and $e(x,t)=E/\rho$. For a nonintegrable system, it is
expected that the evolving system is in local thermal equilibrium and the
three fields $(\rho,v,e)$ should contain the local thermodynamic information
at any space-time point $x,t$. Thus the internal energy per unit mass is
$\epsilon(x,t)=e-{v^{2}}/{2}$, while ideal gas thermodynamics gives
$T(x,t)=2\mu\epsilon(x,t)$ and $P(x,t)=2\rho(x,t)\epsilon(x,t)$.

We consider an initial macrostate where the gas has a finite uniform density
$\rho_{\infty}$, zero flow velocity $v$, and is at zero temperature everywhere
except in a region of width $\sigma$ centered at $x=0$. This is the region of
the blast and we take a smooth Gaussian profile
$E(x,0)=\frac{E}{\sqrt{2\pi\sigma^{2}}}e^{-x^{2}/{2\sigma^{2}}}$. The
procedure to realize this macrostate in the microscopic simulations of the AHP
gas and other details of our numerics are given in the Supplemental Material
supplemental .

(a)(d)

(b)(e)

(c)(f)

Figure 2: (a),(b),(c) Molecular dynamics results for the time evolution of
density, velocity and temperature fields, starting from the initial conditions
corresponding to a Gaussian initial temperature profile and
$\rho(x,0)=\rho_{\infty}=1.5,v(x,0)=0$. The simulation parameters were
$N=L=24000$, $E=32,\mu=1.5$ and an average was taken over $10^{4}$ initial
conditions. (d),(e),(f) This shows the $x/t^{2/3}$ scaling of the data. We
observe a very good collapse of the data at the longest times and an agreement
with the TvNS scaling solution (dashed line).

(a)(d)

(b)(e)

(c)(f)

Figure 3: (a),(b),(c) Evolution of $\rho(x,t),v(x,t),T(x,t)$ obtained from a
numerical solution of the NSF equations (5a)–(5c), starting from the same
initial conditions as used in the simulations for Fig. 2. The other parameters
used in the numerics are $D_{1}=1$, $D_{2}=1$ and $L=4000$. (d),(e),(f)
Scaling plot and comparison with the TvNS solution.

Comparison of simulations with the TvNS solution. – We now present the results
of the microscopic simulations for the evolution of the hydrodynamic fields
$\rho(x,t)$, $v(x,t)$, and $T(x,t)$ starting from the blast-wave initial
conditions. In Fig. 1 we show the spatiotemporal evolution of the three fields
(for individual particle trajectories see Ref. supplemental ). We see a sharp
shock front that evolves sub-ballistically. The mass density and flow velocity
of the gas are peaked around the blast front while the temperature has an
additional peak at the center. In Figs. 2(a),(b),(c) we plot the evolution of
the blast wave at long times. In Figs. 2(d),(e),(f) we show the scaled fields
$\widetilde{\rho}=\rho_{\infty}G-\rho_{\infty},\widetilde{v}=t^{1/3}v(x,t)=(2\alpha/3)\xi
V$ and $\widetilde{T}=t^{2/3}T(x,t)=(4\mu\alpha^{2}/27)\xi^{2}Z$ as functions
of the scaling variable $\xi$. We find an excellent collapse of the data
everywhere except near the blast center. In the region where there is a
collapse of data, we find a perfect agreement with the exact TvNS scaling
solution (plotted as black dashed lines). Close to the origin, the TvNS
solution predicts the singular forms arxiv  $G(\xi)\sim\xi^{1/2}$ and
$Z(\xi)\sim\xi^{-5/2}$ which implies that near the origin the density vanishes
as $\rho\sim x^{1/2}$ and the temperature diverges as $T\sim x^{-1/2}$. These
are unphysical and disagree with simulations. The deviations are caused by
dissipation, specifically heat conduction that becomes important near the
origin. Hence, we need to use the NSF equations.

Comparison of simulations with numerical solution of the NSF equations. — We
now compare simulation results and the TvNS solution with those from the full
dissipative hydrodynamic equations. In one dimension, the NSF equations read

$\displaystyle\partial_{t}\rho+\partial_{x}(\rho v)=0,$ (5a)
$\displaystyle\rho(\partial_{t}+v\partial_{x})v+\partial_{x}(\rho
T/\mu)=\partial_{x}(\zeta\partial_{x}v),$ (5b)
$\displaystyle\frac{\rho^{3}}{2\mu}(\partial_{t}+v\partial_{x})\left(\frac{T}{\rho^{2}}\right)=\partial_{x}(\zeta
v\partial_{x}v+\kappa\partial_{x}T),$ (5c)

where $\zeta$ denotes the bulk viscosity and $\kappa$ is the thermal
conductivity of the system. These transport coefficients can depend on the
fields and, based on the Green-Kubo relations, we expect their temperature
dependence to be $\zeta\sim T^{1/2}$ and $\kappa\sim T^{1/2}$. A recent
numerical study Hurtado2016  suggests the density dependence
$\kappa\sim\rho^{1/3}$. In our numerical study we have thus used the forms
$\zeta=D_{1}T^{1/2}$ and $\kappa=D_{2}\rho^{1/3}T^{1/2}$, where $D_{1}$ and
$D_{2}$ are constants.

We solve Eqs. (5) numerically supplemental ; maccormack1982  for the same
initial conditions as considered in the microscopic simulations, namely
$\rho(x,0)=\rho_{\infty},v(x,0)=0$ and $E(x,0)$ given by the Gaussian form
with total energy $E$. The numerical results are plotted in Figs. 3(a),(b),(c)
which show $\rho(x,t)$, $v(x,t)$ and $T(x,t)$ at different times. In Figs.
3(d),(e),(f) we plot the scaled fields
$\widetilde{\rho},\widetilde{v},\widetilde{T}$ and verify the agreement with
the TvNS solution everywhere except in the core. In the core, the NSF solution
does not have any singularities, unlike the TvNS solution. In Fig. 4 we plot
together the long-time microscopic simulation results and the NSF results and
it can be seen that in the core region, the NSF solution is in better
qualitative agreement with simulations as compared to the TvNS solution.

(a)(b)

Figure 4: Comparison of the scaled fields $\widetilde{\rho},\widetilde{T}$
from microscopic simulations (at time $t=80000$), the NSF equations (at time
$t=6400$), and the vNS scaling solution.

We now estimate the size of the core and outline how to derive the
corresponding scaling solution. Using Eq. (5c) we see that the heat conduction
term becomes important at a length scale $X$ such that
$\rho\,\frac{T}{t}\sim\frac{T^{3/2}}{X^{2}}\,\rho^{1/3}$, where we neglect
constant factors. The TvNS solution in the $\xi\to 0$ limit gives
$\rho\sim\sqrt{\xi}\sim|x|^{1/2}/t^{1/3}$ and $T\sim\xi^{-1/2}\sim
t^{1/3}/|x|^{1/2}$, from which we get an estimate

$X\sim t^{\frac{38}{93}}$ (6)

for the size $X$ of the core where heat conduction is important. The outer
solution (TvNS) and the inner (core) solution are comparable at $|x|=X$ and
this allows us to determine the inner core scaling laws. Thus, the temperature
at the center of the explosion is estimated from $T_{0}\sim
X^{-\frac{1}{2}}\,t^{-\frac{1}{3}}$ and Eq. (6). Similarly, the density at the
center of the explosion can be estimated from $\rho_{0}\sim X^{1/2}/t^{1/3}$
and Eq. (6). We thus arrive at

$\rho_{0}\sim t^{-\frac{4}{31}},\qquad T_{0}\sim t^{-\frac{50}{93}}$ (7)

while the velocity in the core scales as $X/t\sim t^{-{55}/{93}}$.

The hydrodynamic fields in the hot core where dissipative effects play the
dominant role should therefore exhibit scaling behaviors in terms of the
scaled spatial coordinate $\eta=x/X$. For large $x,t$, in a region with
$\eta\sim\mathcal{O}(1)$, we expect the self-similar forms

$\rho=t^{-\frac{4}{31}}\widetilde{G}(\eta),\quad
v=t^{-\frac{55}{93}}\widetilde{V}(\eta),\quad
T=t^{-\frac{50}{93}}\widetilde{Z}(\eta).$ (8)

In Fig. 5 we show that the data from the numerical solution of the NSF
equations approximately satisfy this scaling form, though it appears that the
convergence is somewhat slow. For comparison we also plot the simulation data
at the last two times under the same scaling. A more detailed discussion of
the inner-core scaling solution and its “matching” with the outer-core
solution can be found in Ref. arxiv .

(a)(b)

Figure 5: Verification of the scaling form Eqs. (8), for $\rho$ and $T$, in
the core of the blast. The comparison with data both from simulations and from
the NSF solution are shown. The data are the same as in Figs. 2,3. The value
$\eta=10x/t^{38/93}$ corresponds to $\xi\approx 0.2$.

Conclusions. – We have made a detailed comparison of the predictions of Euler
hydrodynamics for a compressible one-dimensional gas with results from
microscopic dynamics of the hard-core gas. Specifically, we have considered
the blast problem, viz. a localized instantaneous release of energy in a cold
gas. We have derived exact results for the front position and the hydrodynamic
variables and found a remarkable agreement with microscopic simulations.
Deviations were seen in a core region whose size follows the scaling law with
an unusual exponent, $X\sim t^{38/93}$. The position of the shock grows
faster, $R\sim t^{2/3}$, so the relative size of the core region decays to
zero in the long time limit. The width vanishes if measured in units of the
TvNS scaling variable $\xi\sim x/t^{2/3}$. Thus we numerically establish that
the Euler equations provide an asymptotically exact description of the
hydrodynamic behavior of the AHP gas. This agreement is better than in higher
dimensions Barbier2015 ; Barbier2015a ; Barbier2016 ; Joy2017 ; Joy2019 ;
Joy2021 . This is a bit counterintuitive as hydrodynamics is a mean-field
theory that is expected to be more precise in higher dimensions. However, in
one dimension we can consider a truly dilute system as the collisions are
inevitable. Thus we know the exact equation of state in the AHP gas, and we
obtained the TvNS scaling functions exactly. In higher dimensions, one relies
on a virial expansion to get the equation of state and the TvNS solution has
to be numerically found. The low-density ideal gas limit would require much
longer simulation times. We have also checked that unlike in the higher-
dimensional cases Joy2019 ; Joy2021  the local equilibrium, a key assumption
in hydrodynamics, is accurately satisfied in our system supplemental .

The deviation at the core was understood as arising from the contribution of
thermal conduction terms in the energy conservation equation, thus requiring a
study of the full Navier-Stokes-Fourier equation for our 1D gas. For an
accurate comparison with the simulation results, we needed the precise form of
the thermal conductivity, $\kappa$, and in particular its dependence on
temperature and density. For the hard particle gas, with the form
$\kappa\sim\rho^{1/3}T^{1/2}$, our analysis of the NSF equations gives an
estimate $|X|\sim t^{38/93}$ for the core size as well as scaling forms for
the fields which we verified in the numerical solution of the NSF equations.
The inclusion of the dissipative terms leads to a qualitative agreement
between the microscopic simulation results and hydrodynamics even in the core
region, in particular, it cures the erroneous prediction of the TvNS solution,
the divergence of the temperature at the center of the explosion. An open and
challenging problem that remains is to obtain a quantitative agreement between
simulations and hydrodynamics in the core region; this would require the
precise form of the thermal conductivity of our 1D gas. Although we focused on
one dimension where we can provide a detailed numerical comparison, our
approach can be extended to higher dimensions where we find a similar growing
core region where thermal conductivity plays a role. The application of
hydrodynamics to quantum systems is also of much recent interest abanov2011 ;
moore2020  and exploring this would be another interesting direction.

We thank Anupam Kundu, R. Rajesh, Sriram Ramaswamy, Samriddhi Sankar Ray and
Vishal Vasan for useful discussions. We acknowledge support of the Department
of Atomic Energy, Government of India, under Project No. RTI4001.

## References

  * [1] R. K. Zeytounian.  Navier-Stokes-Fourier equations: a rational asymptotic modelling point of view.  Springer Science & Business Media, 2012.
  * [2] R. Esposito and M. Pulvirenti.  From particles to fluids.  In S. Friedlander and D. Serre, editors, Handbook of mathematical fluid dynamics, volume 3, chapter 1, pages 1–82. Elsevier North-Holland, Amsterdam, 2004.
  * [3] I. Gallagher.  From Newton to Navier–Stokes, or how to connect fluid mechanics equations from microscopic to macroscopic scales.  Bull. Amer. Math. Soc., 56(1):65–85, 2019.
  * [4] A. N. Gorban.  Hilbert’s sixth problem: the endless road to rigour.  Phil. Trans. R. Soc. A., 376(2118):20170238, 2018.
  * [5] P. Resibois and M. De Leener.  Wiley, 1977.
  * [6] P. Español, J. G. Anero, and I. Zúñiga.  Microscopic derivation of discrete hydrodynamics.  J. Chem. Phys., 131(24):244117, 2009.
  * [7] H Spohn.  Large scale dynamics of interacting particles.  Springer Science & Business Media, 2012.
  * [8] S. Chapman, T. G. Cowling, and D. Burnett.  The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases.  Cambridge university press, 1990.
  * [9] A. de Masi, R. Esposito, and J. L. Lebowitz.  Incompressible Navier-Stokes and Euler limits of the Boltzmann equation.  Commun. Pure Appl. Math., 42(8):1189–1214, 1989.
  * [10] L. Saint-Raymond.  Hydrodynamic limits of the Boltzmann equation.  Springer, Berlin, 2009.
  * [11] S. Olla, S. R. S. Varadhan, and H-T. Yau.  Hydrodynamical limit for a Hamiltonian system with weak noise.  Commun. Math. Phys., 155(3):523–560, 1993.
  * [12] G.I. Taylor.  The formation of a blast wave by a very intense explosion I. Theoretical discussion.  Proc. R. Soc. A, 201(1065):159–174, 1950.
  * [13] G.I. Taylor.  The formation of a blast wave by a very intense explosion. - II. The atomic explosion of 1945.  Proc. R. Soc. A, 201(1065):175–186, 1950.
  * [14] J. Von Neumann.  John Von Neumann: Collected Works. Theory of games, astrophysics, hydrodynamics and meteorology.  Number 6. Pergamon Press, 1963.
  * [15] L. I. Sedov.  Propagation of strong shock waves.  J. Appl. Math. Mech., 10:241–250, 1946.
  * [16] L. I. Sedov.  Similarity and Dimensional Methods in Mechanics.  Academic Press, New York, 1959.
  * [17] L. D. Landau and E. M. Lifshitz.  Fluid Mechanics.  Pergamon Press, New York, 1987.
  * [18] Y. B. Zel’dovich and Y. P. Raizer.  Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena.  Academic Press, New York, 1967.
  * [19] T. Antal, P. L. Krapivsky, and S. Redner.  Exciting hard spheres.  Phys. Rev. E, 78:030301, Sep 2008.
  * [20] Z Jabeen, R. Rajesh, and P. Ray.  Universal scaling dynamics in a perturbed granular gas.  Europhys. Lett., 89(3):34001, Feb 2010.
  * [21] M. Barbier, D. Villamaina, and E. Trizac.  Blast dynamics in a dissipative gas.  Phys. Rev. Lett., 115:214301, Nov 2015.
  * [22] M. Barbier.  Kinetics of blast waves in one-dimensional conservative and dissipative gases.  J. Stat. Mech., 2015(11):P11019, Nov 2015.
  * [23] M. Barbier, D. Villamaina, and E. Trizac.  Microscopic origin of self-similarity in granular blast waves.  Phys. Fluids, 28(8):083302, 2016.
  * [24] J. P. Joy, S. N. Pathak, D. Das, and R. Rajesh.  Shock propagation in locally driven granular systems.  Phys. Rev. E, 96:032908, Sep 2017.
  * [25] J. P. Joy and R. Rajesh.  Shock propagation in the hard sphere gas in two dimensions: comparison between simulations and hydrodynamics.  preprint arXiv:1907.03416, 2019.
  * [26] J. P. Joy, S. N. Pathak, and R. Rajesh.  Shock propagation following an intense explosion: comparison between hydrodynamics and simulations.  J. Stat. Phys., 182(2):34, 2021.
  * [27] P. L. Garrido, P. I. Hurtado, and B. Nadrowski.  Simple one-dimensional model of heat conduction which obeys Fourier’s law.  Phys. Rev. Lett., 86(24):5486, 2001.
  * [28] A Dhar.  Heat conduction in a one-dimensional gas of elastically colliding particles of unequal masses.  Phys. Rev. Lett., 86:3554–3557, Apr 2001.
  * [29] P Grassberger, W Nadler, and L Yang.  Heat conduction and entropy production in a one-dimensional hard-particle gas.  Phys. Rev. Lett., 89:180601, Oct 2002.
  * [30] G. Casati and T. Prosen.  Anomalous heat conduction in a one-dimensional ideal gas.  Phys. Rev. E, 67(1):015203, 2003.
  * [31] P Cipriani, S Denisov, and A Politi.  From anomalous energy diffusion to levy walks and heat conductivity in one-dimensional systems.  Phys. Rev. Lett., 94(24):244301, 2005.
  * [32] S. Chen, J. Wang, G. Casati, and G. Benenti.  Nonintegrability and the Fourier heat conduction law.  Phys. Rev. E, 90(3):032134, 2014.
  * [33] P. I. Hurtado and P. L. Garrido.  A violation of universality in anomalous Fourier’s law.  Sci. Rep., 6(1):38823, 2016.
  * [34] H. Zhao and W. Wang.  Fourier heat conduction as a strong kinetic effect in one-dimensional hard-core gases.  Phys. Rev. E, 97:010103, Jan 2018.
  * [35] S. Lepri, R. Livi, and A. Politi.  Too close to integrable: Crossover from normal to anomalous heat diffusion.  Phys. Rev. Lett., 125:040604, Jul 2020.
  * [36] Y. Du, H. Li, and L. P. Kadanoff.  Breakdown of hydrodynamics in a one-dimensional system of inelastic particles.  Phys. Rev. Lett., 74:1268–1271, Feb 1995.
  * [37] P. I. Hurtado.  Breakdown of hydrodynamics in a simple one-dimensional fluid.  Phys. Rev. Lett., 96(1):010601, 2006.
  * [38] C. B. Mendl and H. Spohn.  Shocks, rarefaction waves, and current fluctuations for anharmonic chains.  J. Stat. Phys., 166(3-4):841–875, Oct 2017.
  * [39] X. Cao, V. B. Bulchandani, and J. E. Moore.  Incomplete thermalization from trap-induced integrability breaking: Lessons from classical hard rods.  Phys. Rev. Lett., 120:164101, Apr 2018.
  * [40] G Casati and T Prosen.  Mixing property of triangular billiards.  Phys. Rev. Lett., 83(23):4729, 1999.
  * [41] See supplemental material for technical details.
  * [42] G. I. Barenblatt.  Scaling, Self-similarity, and Intermediate Asymptotics: Dimensional Analysis and Intermediate Asymptotics.  Cambridge Texts in Applied Mathematics. Cambridge University Press, 1996\.
  * [43] S. Ganapa, S. Chakraborti, P. L. Krapivsky, and A. Dhar.  The Taylor - von Neumann - Sedov blast-wave solution: comparisons with microscopic simulations of a one dimensional gas.  preprint arXiv:2010.15868, 2020.
  * [44] R. W. MacCormack.  A numerical method for solving the equations of compressible viscous flow.  AIAA journal, 20(9):1275–1281, 1982.
  * [45] J. A. Joseph, J. E. Thomas, M. Kulkarni, and A. G. Abanov.  Observation of shock waves in a strongly interacting fermi gas.  Phys. Rev. Lett., 106(15):150401, 2011.
  * [46] V. B. Bulchandani, C. Karrasch, and J. E. Moore.  Superdiffusive transport of energy in one-dimensional metals.  PNAS, 117(23):12713–12718, 2020.

## Supplemental material for ‘Blast in a One-Dimensional Cold Gas: From
Newtonian Dynamics to Hydrodynamics’

##  .1 Exact TvNS scaling solution of Euler equations

We start with the equations

$\displaystyle\partial_{t}\rho+\partial_{x}(\rho v)=0,$ (9a)
$\displaystyle\partial_{t}(\rho v)+\partial_{x}(\rho v^{2}+P)=0,$ (9b)
$\displaystyle\partial_{t}(\rho e)+\partial_{x}(\rho ve+Pv)=0,$ (9c)

where we recall that $e=v^{2}/2+T/(2\mu)$ and $P=\rho T/\mu$. We seek a self-
similar scaling solution of these equations with a shock at the location
$R(t)$. At the shock the condition of weak solutions lead to the the following
Rankine-Hugoniot boundary conditions

$\displaystyle\frac{\rho(R)v(R)}{\rho(R)-\rho_{\infty}}=U,$ (10)
$\displaystyle\frac{\rho(R)v^{2}(R)+P(R)}{\rho(R)v(R)}=U,$ (11)
$\displaystyle\frac{\rho(R)v(R)e(R)+P(R)v(R)}{\rho(R)e(R)}=U,$ (12)

where $U=\dot{R}$. In addition, we need to use the condition of energy
conservation:

$\displaystyle\int_{0}^{R}dx\rho(v^{2}+T/\mu)=E.$ (13)

We seek the hydrodynamic variables in the scaling form

$\rho=\rho_{\infty}G(\xi),\quad v=\frac{2x}{3t}\,V(\xi),\quad
T=\frac{4\mu}{27}\,\frac{x^{2}}{t^{2}}\,Z(\xi),$ (14)

where the re-scaled fields depend on the single dimensionless variable

$\xi=\frac{x}{R}\,,\quad{\rm
where}~{}R(t)=\left(\frac{Et^{2}}{A\rho_{\infty}}\right)^{1/3}.$ (15)

Using the scaling form, the value $\gamma=3$ of the adiabatic index for the
one-dimensional gas of point particles, and after some simplification, the
Rankine-Hugoniot boundary conditions lead to

$G(1)=2,\quad V(1)=\tfrac{1}{2},\quad Z(1)=\tfrac{3}{4}.$ (16)

Plugging the scaling forms (14)–(15) into Eqs. (9a)–(9b) and defining
$\ell=\ln\xi$, we obtain

$\displaystyle\frac{dV}{d\ell}+(V-1)\,\frac{d\ln G}{d\ell}=-V,$ (17)
$\displaystyle\frac{d\ln Z}{d\ell}-2\,\frac{d\ln G}{d\ell}=\frac{3-2V}{V-1}.$
(18)

Using the fact that energy is conserved in the region $(-\xi R(t)<x<\xi
R(t))$, one can express [17, 43] $Z$ through $V$

$Z=\frac{3(1-V)V^{2}}{(3V-1)}.$ (19)

The total energy conservation condition in Eq. (13) gives

$\displaystyle
E=\rho_{\infty}\,\frac{4}{9}\,\frac{R^{3}}{t^{2}}\int_{0}^{1}d\xi\,\xi^{2}\,\frac{2V^{3}}{(3V-1)}\,G,$
(20)

where we have used (19). Combining (20) with (15) we obtain

$A=\frac{4}{9}\,\int_{0}^{1}d\xi\,\xi^{2}\,\frac{2V^{3}}{3V-1}\,G.$ (21)

As shown in [43] it is possible to obtain an explicit solution of the set of
equations (17,18,19) with the boundary condition (16):

$\displaystyle\xi^{5}=2^{-\frac{4}{3}}\,(3V-1)^{2}\,V^{-\frac{10}{3}}\,(3-4V)^{-\frac{11}{3}}\,,$
(22a) $\displaystyle
G=2^{\frac{16}{5}}\,(1-V)^{2}\,(3V-1)^{\frac{1}{5}}\,(3-4V)^{-\frac{11}{5}}.$
(22b)

Plugging this solution into Eq. (21) we reduce the integral over $\xi$ to an
integral over $V$ which is computed to give $A=152/1071$. Thus we have a
completely closed form expression for the TvNS solution for the 1D ideal gas.

##  .2 Specification of microscopic dynamics and initial conditions

Let $q_{j}$, $p_{j}$, and $m_{j}$ be respectively the position, momentum, and
mass of the $j^{th}$ particle with $j=1,2,\ldots,N$. For our binary mass
model, we choose all even numbered particles to have mass $m_{1}$ and all odd
ones $m_{2}$. Between collisions the particles move ballistically with
constant speeds. After collision between particles $i$ and $i+1$ moving with
velocities $v_{i}$ and $v_{i+1}$, their post-collisional velocities follow
from the conservation of energy and momentum:

$\displaystyle v_{i}^{\prime}$
$\displaystyle=\frac{[2m_{i+1}v_{i+1}+v_{i}(m_{i}-m_{i+1})]}{(m_{i}+m_{i+1})}\,,$
(23) $\displaystyle v_{i+1}^{\prime}$
$\displaystyle=\frac{[2m_{i}v_{i}+v_{i+1}(m_{i+1}-m_{i})]}{(m_{i}+m_{i+1})}\,.$
(24)

In our numerical simulations we consider the following Gaussian profile for
the initial energy:

$\displaystyle
E(x,0)=\frac{E}{\sqrt{2\pi\sigma^{2}}}e^{-x^{2}/{2\sigma^{2}}}.~{}$ (25)

The total energy of the blast is $E$.

To obtain this in the particle model, we first distribute $N$ particles
labelled $i=-N/2+1,-N/2+2,\ldots,N/2$ uniformly between $x=-L/2$ to $L/2$, so
that the number density is $n_{\infty}=N/L=\rho_{\infty}/\mu$, where
$\mu=(m_{1}+m_{2})/2$ is the mean mass and $\rho_{\infty}$ the ambient mass
density of the gas. For the $N_{c}$ particles around the center of the
explosion with labels $i=-N_{c}/2+1,-N_{c}/2+2,\ldots N_{c}/2$, we choose
their velocities from the Maxwell distribution with equal temperature:
$\text{Prob}(v_{i})=\sqrt{m_{i}/(2\pi T_{c})}e^{-m_{i}v_{i}^{2}/(2T_{c})}$
with $T_{c}=2E/N_{c}$. The velocities of other particles are set to zero. The
size of the initial blast is thus approximately $s=N_{c}/n_{\infty}$.

We note that
$E(x,0)=\langle\sum_{i=1}^{N}m_{i}v_{i}(0)^{2}/2\delta(x-q_{i}(0))\rangle$.
Since, for large $N$, each particle’s position is a Gaussian with mean
$\bar{q}_{i}=i/n_{\infty}$ and variance $\sigma^{2}=L/(4n_{\infty})$, we then
get:

$\displaystyle E(x,0)$
$\displaystyle=\sum_{i=-N_{c}/2+1}^{N_{c}/2}\frac{T_{i}}{2\sqrt{2\pi\sigma^{2}}}e^{-(x-\bar{q}_{i})^{2}/(2\sigma^{2})}$
$\displaystyle\approx\frac{n_{\infty}T_{i}}{2}\int_{-s/2}^{s/2}dy\frac{e^{-(x-y)^{2}/(2\sigma^{2})}}{\sqrt{2\pi\sigma^{2}}}$
$\displaystyle=\frac{n_{\infty}T_{i}}{4}\left[\operatorname{erf}{\left(\frac{s-2x}{2\sqrt{2}\sigma}\right)}+\operatorname{erf}{\left(\frac{s+2x}{2\sqrt{2}\sigma}\right)}\right].$
(26)

When $s\ll\sigma$ this simplifies to

$\displaystyle E(x,0)\approx
E\,\frac{e^{-x^{2}/(2\sigma^{2})}}{\sqrt{2\pi\sigma^{2}}}.$ (27)

We have verified that our long-time results are independent of the initial
profile as long as it is localized and the total energy $E$ is fixed [43].

Figure 6: This plot shows the trajectories of $N=1000$ particles, starting
from initial conditions corresponding to a Gaussian initial temperature
profile and $\rho(x,0)=\rho_{\infty}=1.5,v(x,0)=0$. The expected shock front
position $R(t)$ is shown as a dashed line. The stationary particles are not
shown.

The hydrodynamic fields are obtained from the microscopic variables using the
standard relations:

$\displaystyle\rho(x,t)$
$\displaystyle=\sum_{j=1}^{N}m_{j}\left\langle\delta(q_{j}(t)-x)\right\rangle,$
(28a) $\displaystyle p(x,t)$
$\displaystyle=\rho(x,t)v(x,t)=\sum_{j=1}^{N}m_{j}\left\langle
v_{j}\delta(q_{j}(t)-x)\right\rangle,$ (28b) $\displaystyle E(x,t)$
$\displaystyle=\rho(x,t)e(x,t)=\sum_{j=1}^{N}\frac{m_{j}}{2}\left\langle
v_{j}^{2}\delta(q_{j}(t)-x)\right\rangle.$ (28c)

Here $\langle...\rangle$ indicates an average over an initial distribution of
microstates that correspond to the same initial macrostate.

##  .3 Details of numerical techniques

Molecular dynamics simulations: In all our simulations we took
$m_{1}=1,m_{2}=2$ (so that $\mu=(m_{1}+m_{2})/2=1.5$), $\rho_{\infty}=1.5$ and
$E=32,N_{c}=32$. In our largest simulations we took $N=24000$, $L=24000$ and
averaged over an ensemble of ${\cal{R}}=10^{4}$ initial conditions. For each
microscopic initial condition, we evolved the system with the Hamiltonian
dynamics. The molecular dynamics for this system was done using an event-
driven algorithm that updates the system between successive collisions.

Solution of the NSF equations: Our numerical solution relies on the MacCormack
method [44], which is second order in both space and time. We have used a
discretization $dx=0.1$ and $dt=0.001$. The diffusion constants in the
dissipative terms were set to the values $D_{1}=D_{2}=1$ and the system size
was taken as $L=4000$. We have evolved the system up to time $t=6400$ which is
before the energy reaches the boundary. We have checked numerically the
convergence of the method by using different values of $dx$ and $dt$. Note
that for the hyperbolic Euler equations, with no dissipative terms, the
numerical approach leads to strong dispersive effects near the shock front.
The presence of $D_{1}$ and $D_{2}$ terms in the equations, in addition to
being necessary to explain the microscopic observations, also ensure that the
solution does not develop dispersion and is stable at long times.

Figure 7: The spatial profile of the skewness defined via
$S=n^{-1}\sum_{i}[v_{i}-u(x,t)]^{3}\delta(x_{i}-x)$ and the kurtosis defined
via $\chi=n^{-1}\sum_{i}m_{i}^{2}[v_{i}-u(x,t)]^{4}\delta(x_{i}-x)$ are
plotted versus the scaled spatial coordinate $\xi=x/R(t)$. Here $n$ is the
local particle number density. The skewness is scaled by $[{\rm
thermal~{}velocity}]^{3}$ while the curtosis is scaled by the expected thermal
value of $3T^{2}(x,t)$. The results at two different times are shown.

##  .4 Probing the validity of the hydrodynamic description

Hydrodynamics is expected to provide an adequate description for phenomena
where the characteristic macro scale ($R$ in our case) greatly exceeds the
micro scale (separations between adjacent particles in 1D). It is thus not
expected to describe the structure of the shock wave since the width of the
shock is its scale and the width is of the order of the separation (the micro
scale).

The local equilibrium (LE) assumption is an important requirement for the
validity of the hydrodynamic description and verifying this would be a direct
probe of the validity of the hydrodynamic description. Previous studies in two
and three dimensions [25, 26] find strong deviations from the LE — the
observed skewness and kurtosis of the local velocity distributions indicate
non-Gaussianity. The lack of the LE can be a reason for the observed
disagreement with the hydrodynamic (TvNS) predictions. Here we show results
for these quantities obtained from MD simulations of our 1D system. Fig. 7
shows the LE assumption is satisfied quite accurately, significantly better
than in higher dimensions [25, 26]. We do see small deviations from the LE
which, interestingly, peak in the core and shock regions. While this could
provide a possible explanation for the disagreement between the NSF and
simulations in the core region, we note however that we continue to have good
agreement between these theories near the shock front.

However, note that the structure of the shock wave is certainly not accounted
by Euler hydrodynamics. Studying this and also the approach to the predictions
of hydrodynamics is an avenue for future work. For example, finding the more
precise formula for the position of the shock wave is an open problem.
Hydrodynamics tells us that

$\displaystyle R=[1071Et^{2}/(152\rho)]^{1/3}$ (29)

but the deterministic hydrodynamic equations say nothing about the sub-leading
term. The role of higher order terms in the Chapman-Enskog-Burnett gradient
expansion [8], towards determining corrections to the leading scaling form,
also need to be probed further.