# Connect the Dots: In Situ 4D Seismic Monitoring of CO2 Storage with Spatio-
temporal CNNs

Shihang Feng2,∗, Xitong Zhang1,3, Brendt Wohlberg2, Neill Symons1, and Youzuo
Lin1,∗ 1: Earth and Environmental Sciences Division, Los Alamos National
Laboratory, Los Alamos, NM, 87544 USA.2: Theoretical Division, Los Alamos
National Laboratory, Los Alamos, NM, 87544 USA.3: Department of Computational
Mathematics, Science and Engineering, Michigan State University, East Lansing,
MI 48824 USA.∗Correspondence to: S. Feng<EMAIL_ADDRESS>and Y. Lin
(ylin@lanl.gov).

###### Abstract

4D seismic imaging has been widely used in CO2 sequestration projects to
monitor the fluid flow in the volumetric subsurface region that is not sampled
by wells. Ideally, real-time monitoring and near-future forecasting would
provide site operators with great insights to understand the dynamics of the
subsurface reservoir and assess any potential risks. However, due to obstacles
such as high deployment cost, availability of acquisition equipment, exclusion
zones around surface structures, only very sparse seismic imaging data can be
obtained during monitoring. That leads to an unavoidable and growing knowledge
gap over time. The operator needs to understand the fluid flow throughout the
project lifetime and the seismic data are only available at a limited number
of times. This is insufficient for understanding the reservoir behavior. To
overcome those challenges, we have developed spatio-temporal neural-network-
based models that can produce high-fidelity interpolated or extrapolated
images effectively and efficiently. Specifically, our models are built on an
autoencoder, and incorporate the long short-term memory (LSTM) structure with
a new loss function regularized by optical flow. We validate the performance
of our models using real 4D post-stack seismic imaging data acquired at the
Sleipner CO2 sequestration field. We employ two different strategies in
evaluating our models. Numerically, we compare our models with different
baseline approaches using classic pixel-based metrics. We also conduct a blind
survey and collect a total of 20 responses from domain experts to evaluate the
quality of data generated by our models. Via both numerical and expert
evaluation, we conclude that our models can produce high-quality 2D/3D seismic
imaging data at a reasonable cost, offering the possibility of real-time
monitoring or even near-future forecasting of the CO2 storage reservoir.

###### Index Terms:

4D Seismic Imaging, In Situ Monitoring, CO2 Sequestration, Data Interpolation
& Extrapolation, Data Scarcity, Spatial-temporal Analysis

##  I Introduction

To assess the safety and effectiveness of CO2 sequestration projects,
effective monitoring techniques are required to build a 4D (3D in space and 1D
in time) digital subsurface models that accurately matches the CO2 injection
over a significant period of time (e.g. 10 - 200 years). Among all available
geophysical monitoring techniques (such as seismic, EM, gravity, ERT, etc.),
4D seismic imaging is one of the most effective and reliable techniques for
monitoring and verification of CO2 sequestration projects [1]. However, due to
insufficient data, in-situ 4D seismic monitoring of the reservoir changes is
technically challenging or even infeasible.

The data scarcity is mainly caused by the practical limitations of high cost,
availability of equipment, exclusion zones around surface structures, etc.,
that makes it difficult to deploy sufficient frequency of data acquisition
[2]. For example, consider the world’s first industrial-scale CO2 storage
operation project, Sleipner, a well-studied area in the North Sea [3]. Since
October 1996, Equinor, the site operator, has injected about 1 million tonnes
(Mt) of super-critical CO2 per year [4]. In order to monitor the migration of
injected CO2 and plume geometry development, comprehensive 3D-time lapse
seismic surveys have been carried out. Due to practical limitations, data were
acquired only in a limited number of years (in 1994, 1999, 2001, 2002, 2004,
2006, 2008, 2010, 2013, and 2016). The missing data in between creates an
unavoidable knowledge gap, which hinders the in-situ monitoring of the CO2
migration. To improve the security and applicability of the CO2 sequestration
project, these missing data are required to be estimated for real-time
monitoring, which is technically challenging.

Recently, similar challenges and tasks have been actively studied in the
computer vision community for generating new samples from a sequence of
temporally-dependent images, namely, “in-between image interpolation” [5, 6,
7]. Broadly speaking, there are two categories of interpolation approaches,
i.e., those based on pixels and those based on latent representations. Linear
interpolation is a typical example of a pixel-based approach. Pixel-based
methods usually cause problems when input images are temporally far apart,
which would lead to a loss of temporal dependence and create errors like
cross-dissolve effects [8]. On the other hand, by leveraging semantic
information, the latent representation-based approaches alleviate the problems
of pixel-based models to some extent. However, directly applying those
existing in-between image interpolation techniques to our 4D seismic imaging
data is not trivial due to the differences in the nature of the problems as
well as the spatio-temporal characteristics.

4D post-stack seismic imaging (i.e., time-lapse 3D imaging) leverages
repeating seismic surveys in succession. Via time-lapse 3D imaging techniques,
time shifts and amplitude changes in the data can be measured and further
interpreted as the changes in the reservoir and overburden [9, 10, 11]. In a
typical CO2 sequestration project such as Sleipner, the seismic data can
change rapidly as the CO2 injection proceeds. This is due to the large
perturbation of subsurface velocity as the fluid changes from pure brine to
brine with low CO2 saturation, while the velocity changes much less as CO2
saturation increases further [12]. 4D post-stack seismic imaging data can
simultaneously detect spatial and temporal changes [13] if the sampling is
sufficiently dense. Spatially, the CO2 plume in post-stack imaging with the
Sleipner data is characterized as several bright sub-horizontal seismic
reflections within the reservoir unit. The seismic reflections mostly come
from thin layers of CO2 trapped beneath very thin intra-reservoir mudstones
and the reservoir caprock and the push-down is caused by the decrease of
seismic velocity due to CO2-saturated rock. Temporally, the migration and
development of CO2 plume over time is characterized by repeating surveys at
relatively sparse times. Hence, filling the knowledge gap for post-stack
seismic imaging would need to respect the spatio-temporal characteristics of
the existing sparse measurements without violating the physics. In the problem
in-situ monitoring, we need to generate in-between images (interpolation) and
also to forecast new images from a few past measurements (extrapolation) which
is more technically challenging.

To overcome these aforementioned challenges, we propose a deep neural network-
based interpolation/extrapolation model incorporating domain knowledge to
“connect the dots” between the seismic data in years. Our method builds on an
autoencoder to obtain latent variables of input sequential images. To further
enhance the temporal correlation, we not only incorporate the long short-term
memory (LSTM) [14] but also use a new loss function regularized with optical
flow [15]. Particularly, the LSTM structure is used to recognize the
sequential and temporal dynamics, which are leveraged for forecasting. The
optical flow describes the spatial shifting of objects in time-varying images.
We use it as a regularization term to provide an additional layer of
constraints when generating new images.

We validate the performance of our techniques using both real 2D and 3D post-
stack seismic imaging data acquired from the Sleipner project. Similar to
previous work, we first compare our approaches to baseline methods by using
standard per-pixel error metrics. However, classic per-pixel measures are
insufficient for assessing structured data such as images [16]. We carry out
an expert evaluation test to further justify the quality of our generated
images. Via both standard metrics and expert evaluation, we show that our
model outperforms existing approaches and effectively and efficiently produces
realistic images that could facilitate in-situ monitoring of CO2 changes.

This paper is organized into the following sections. After the introduction in
the first section, we provide related work in data interpolation from both the
seismic exploration and computer vision communities. We then present the
foundation and theory including 3D time-lapse seismic imaging and Sleipner
data set, autoencoder, and optical flow in the third section. Our new proposed
interpolation and extrapolation networks are provided in the fourth section.
The fifth section shows numerical results using the Sleipner data. A
discussion of our interpolation and extrapolation networks is presented in the
sixth section. Finally, the conclusion is given in the last section.

##  II Related Works

###  II-A Interpolation and Extrapolation in Seismic Exploration

Seismic interpolation is an important research topic in the field of seismic
data processing. The prediction filters in the F-X domain [17], support vector
regression (SVR) [18, 19], U-net [20] and CycleGAN [21] had been applied to
reconstruct seismic data from under-sampled or missing traces in the trace
interpolation. Interpolation methods have been studied to obtain continuous
information between 4D seismic surveys. Linear and Lagrange interpolation of
the seismic response were used for calculating the onset time map [22]. An
amplitude interpolation method based on a probability formulation was
developed to estimate the seismic data between the infrequent seismic survey
[23]. Although these methods can produce smooth transitions between the
repeated survey, interpolation uncertainty increases with the time gap between
consecutive seismic surveys.

Prediction of seismic data is also important in the study of 4D seismic
survey. A comprehensive analysis of the geophysical data and geology is
usually required to increase the confidence of the prediction [24, 25]. The
amplitude interpolation method had been extended to predict the seismic data,
but a large uncertainty exists with the extrapolation [23]. With some physical
regularization such as optical flow, we can mitigate the uncertainty problem
in our extrapolation network.

TABLE I: The total amount of CO2 injected volume along with the time | Year
---|---
% | 1994 | 1999 | 2001 | 2002 | 2004 | 2006 | 2008 | 2010
Mass (Mt) | 0 | 2.3 | 4.2 | 5.1 | 6.9 | 8.4 | 10.1 | 12.0

###  II-B Interpolation in Computer Vision

Image interpolation has been the subject of intense study in computer vision.
Existing methods can be generally categorized as pixel-based or latent-
representation based [5]. The former performs a simple manipulation of pixels
by considering the motion of the objects from the pixel level. The latter
learns a feature-based representation of object motion from the semantic level
by convolving input images with spatially adaptive kernels.

Pixel-based interpolation between images introduces ghosting artifacts, a sign
of departure from the underlying manifold [26]. Two directions had been
studied to overcome the limitations of those approaches. The first direction
is the motion-compensated frame methods [27, 28] and the second is the phase-
based methods [29, 30]. The motion-compensated approaches are based on the
assumption that motion in images is smooth and continuous. They estimate the
motion based on previous and current images so that the in-between images can
be generated by averaging the pixels in the images indicated by motion vectors
[27, 28]. The phase-based methods assume that small motion can be encoded
using phase shift [29, 30].

With recent advances in neural networks, the latent-representation-based
approaches had shown promising results in the task of in-between image
interpolation [31, 32, 33, 34, 35]. Due to the capability of disentangling
important factors of variation in the dataset, those latent-representation-
based approaches are particularly effective in producing semantically smooth
interpolations [35, 33]. A critic was introduced in the training of the
interpolation network to generate more realistic new samples [32]. Deep neural
networks and domain knowledge were combined to improve the performance of the
interpolation [34].

##  III Background

###  III-A 4D Seismic Imaging and Sleipner Geology

4D post-stack seismic imaging for monitoring CO2 sequestration involves
repeated seismic surveys in time-lapse to image the subsurface changes due to
CO2 injection [36] over time. It provides an understanding of the spatial
movement and distribution of hydrocarbon and fluids, including supercritial
CO2 and the saline water saturated by CO2, over time. The seismic data with an
initial baseline survey before the injection and several follow-up surveys
after a few years of injection are collected and processed with the same
processing parameters. The anomalies between the processed data with the
baseline survey and the follow-up surveys are attributable to changes in the
reservoir [37], such as fluid movement [38] and CO2 injection [39].

The subsurface geology of the Sleipner field in the North Sea has been well
studied and shown to be ideal for CO2 sequestration [40]. Particularly, the
storage reservoir is a saline aquifer of regional extent stretching for more
than 400 km north to south and 50$\sim$100 km east to west, locally exceeding
thicknesses of 300 m [41] shown in Figure 1. CO2 is injected into highly
porous Utsira Formation, a saline aquifer approximately 1 km below the surface
and 200 km from the coast [42]. The injection is conducted through a single
near-horizontal injection well at a depth of 1,012 m below sea level (bsl)
[43]. 4D time-lapse seismic surveys were implemented from 1994 to 2016. One
pre-injection dataset and ten post-injection datasets are collected in a 4 km
$\times$ 7 km area. The collected data are processed to generate 3D seismic
stacked sections [44]. With the injection of CO2, a growing CO2 plume is
manifested as a series of bright reflections in the sections [40]. The
reflections mostly arise from thin layers of CO2 trapped beneath very thin
intra-reservoir mudstone and the reservoir caprock [45]. Moreover, the plume
also causes a velocity push-down since the seismic waves travel much slower
through CO2-saturated rock than through the virgin aquifer [4].

Figure 1:  a) Location map [12]. b) Utsira two-way traveltime (TWT) Map. Red
is shallow, blue is deep [46].

Repeated 3D seismic data were collected in 1994, 1999, 2001, 2002, 2004, 2006,
2008, 2010, 2013, and 2016 for monitoring. However, only data from 1999 to
2010 are openly released. The total volume of injected CO2 in each year is
shown in Table LABEL:tab:inject. The dataset in each year is a 3D post-stack
seismic section with 2 seconds recording in a time spacing of 0.002 s. Each
section has 480 grids in the inline direction and 144 grids in the crossline
direction. The data differences caused by the CO2 injection are mostly after
0.5 s. We use a Tesla V100 GPU in the training of the networks. Using the full
temporal length (i.e., from 0.5 s to 2.0 s) will, unfortunately, exceed the
GPU memory limit. Via numerous trying and testing, we found out that an image
data with a temporal length from 0.5 s to 1.56 s is the maximum allowable size
that not only fits right in our GPU device, but also contains all the critical
information (i.e., the overburden, the new reflections generated by CO2 plume,
and the shift of the reflection events). As a result, we extract the seismic
sections from 0.5 to 1.56 s, resulted in the size of the extracted seismic
sections as 480$\times$480$\times$144\. A total of 8 extracted seismic
sections from 1994 to 2010 are available for this study.

###  III-B Autoencoders and Latent Representation

The autoencoder is one of the prevalent unsupervized methods to generate a
compressed representations of its input [47]. The representation is usually
exploited to discover the underlying structural patterns of training data in
the low-dimensional latent space by downstream machine learning models. In
general, the autoencoder consists of two parts: an encoder network and a
decoder network. The encoder network aims to map the data from high-
dimensional input to a lower-dimensional latent space. The decoder network
recovers the data from the latent space. The encoder and decoder networks are
trained simultaneously to minimize reconstruction errors. Unlike conventional
dimension reduction methods, e.g., principal component analysis (PCA), it is a
non-linear compression model transforming inputs to a lower-dimensional latent
space, which can capture the complex distribution of training data. The
compressed representation can be reconstructed back to the original data space
that ensures that the representation includes the critical information of
inputs.

Figure 2: A schematic illustration of the interpolation network which contains
two encoder networks and three decoder networks. The observed data at time $n$
and $n+m$ are fed into the encoder networks to extract the latent spaces. Then
the latent spaces and their linear combination are fed in the decoder networks
to generate the data at $n$, $n+i$ and $n+m$.

###  III-C Optical Flow

Optical flow is a technique, originating in computer vision, to detect the
motion of objects between consecutive frames of the sequence [15, 48]. It
consists of vectors at each pixel location that represent the motions of the
pixel from one frame to another. Particularly, given two successive images
$P(x,y,z,t_{1})$ and $P(x,y,z,t_{2})$, we need to find the motion vector
$\mathbf{m}(x,y,z)$ from $P(x,y,z,t_{1})$ and $P(x,y,z,t_{2})$. In 4D seismic
imaging problem, $P(x,y,z,t_{1})$ and $P(x,y,z,t_{2})$ are seismic data at two
close sequential times. A point at location $(x,y,z)$ in $P(x,y,z)$ will move
by $\Delta{x}$, $\Delta{y}$ and $\Delta{z}$ after the time interval
$\Delta{t}$ and this can be represented as:

$P(x,y,z,t)=P(x+\Delta{x},y+\Delta{y},z+\Delta{z},t+\Delta{t}).$ (1)

Assuming the movement to be small, the image at
$P(x+\Delta{x},y+\Delta{y},z+\Delta{z},t+\Delta{t})$ can be expanded using a
Taylor series:

$\displaystyle P(x+\Delta{x},y+\Delta{y},z+\Delta{z},t+\Delta{t})=$
$\displaystyle
P(x,y,z,t)+P_{x}\Delta{x}+P_{y}\Delta{y}+P_{z}\Delta{z}+P_{t}\Delta{t}+...,$
(2)

where $P_{x}$, $P_{y}$, $P_{z}$ and $P_{t}$ are partial derivative with
respect to $x$, $y$, $z$ and $t$. From Eqs. (1) and (2), ignoring the higher
order terms (i.e.,“…” in Eq. (2)), it follows that

$P_{x}u+P_{y}v+P_{z}w+P_{t}=0,$ (3)

where $u$, $v$ and $w$ are the $x$, $y$, $z$ components of the motion vector
$\mathbf{m}(x,y,z)$, which are defined as

$u=\frac{\Delta{x}}{\Delta{t}},v=\frac{\Delta{y}}{\Delta{t}},w=\frac{\Delta{z}}{\Delta{t}}.$
(4)

Combining the $\ell_{2}$ loss function with smoothing regularization, the
objective function can be written as

$\zeta^{\mathrm{op}}=\iiint(P_{x}u+P_{y}v+P_{z}w+P_{t})^{2}+\lambda^{2}\,M^{2}\,dxdydz,$
(5)

where $\lambda$ is the weight for the smoothing regularization, and $M$ is the
magnitude of the flow gradient given by

$M({u},{v},{w})=\sqrt{\|\nabla{u}\|^{2}+\|\nabla{v}\|^{2}+\|\nabla{w}\|^{2}}.$
(6)

An iterative solution can be obtained as

$u^{\mathrm{k+1}}=\bar{u}^{\mathrm{k}}-\frac{P_{x}({P_{x}\bar{u}^{\mathrm{k}}+P_{y}\bar{v}^{\mathrm{k}}+P_{z}\bar{w}^{\mathrm{k}}+P_{t}})}{\alpha_{op}^{2}+P_{x}^{2}+P_{y}^{2}+P_{z}^{2}},$
(7)
$v^{\mathrm{k+1}}=\bar{v}^{\mathrm{k}}-\frac{P_{y}({P_{x}\bar{u}^{\mathrm{k}}+P_{y}\bar{v}^{\mathrm{k}}+P_{z}\bar{w}^{\mathrm{k}}+P_{t}})}{\alpha_{op}^{2}+P_{x}^{2}+P_{y}^{2}+P_{z}^{2}},$
(8)
$w^{\mathrm{k+1}}=\bar{w}^{\mathrm{k}}-\frac{P_{z}({P_{x}\bar{u}^{\mathrm{k}}+P_{y}\bar{v}^{\mathrm{k}}+P_{z}\bar{w}^{\mathrm{k}}+P_{t}})}{\alpha_{op}^{2}+P_{x}^{2}+P_{y}^{2}+P_{z}^{2}},$
(9)

where $\bar{u}$, $\bar{v}$ and $\bar{w}$ are the local spatially smoothed
values of $u$, $v$ and $w$, respectively. $k$ and $k+1$ are the iteration
numbers and $\alpha_{op}$ is a hyper-parameter that helps avoid division by
zero.

##  IV Methods

Figure 3: A schematic illustration of the extrapolation network which contains
four encoder networks and five decoder networks. The observed data change from
time $n$ to $n+3i$ are fed into the encoder networks to extract the latent
spaces. The latent spaces are connected using an LSTM network to predict the
latent space at time $n+4i$. All these latent spaces are fed in the decoder
networks to predict the data change from $n$ to $n+4i$.

###  IV-A Method Overview

In this section, we develop two different networks for interpolation and
extrapolation. For the interpolation task, inspired by the variational
autoencoder [49], we assume a linear pattern among samples in the latent space
when the CO2 volume increases over time. Thus, time sequences of samples can
be generated by linearly interpolating features in the latent space. The
extrapolation task is much more challenging in the sense that the problem is
more severely ill-posed than interpolation. Prior knowledge needs to be
exploited to make the problem well-posed so that reasonable solutions may be
obtained. Here, we consider the consistency of the physics of the CO2
migration throughout the monitoring, which leads us to incorporate the LSTM
structure and the optical-flow regularization with our extrapolation network
model.

###  IV-B Interpolation Network

Due to the lack of sufficient observations of the underground CO2 plume for
training our interpolation network, it is desirable to introduce
regularization during training to alleviate the overfitting risk. One
reasonable regularization is to assume that the representations lie on the
same line in the latent space. Although the migration patterns are actually
nonlinear, our interpolation network can project time-series observations onto
a line by nonlinear transformation with multiple convolutional layers.

Our interpolation network consists of two encoder networks and three decoder
networks, as shown in Figure 2. Both the encoders and the decoders share their
weights with other encoders and decoders. The observed seismic data with total
CO2 injected volume $n$ and $n+m$, $d_{n}$ and $d_{n+m}$, are fed to the
encoders to obtain their latent space representations $L_{n}$ and $L_{n+m}$.
Based on the linear regularization in the latent space, the interpolated
latent space $L$ with total CO2 injected volume $n+i$ can be represented as

$L_{n+i}=L_{n}+\frac{m-i}{m-n}(L_{n+m}-L_{n}),$ (10)

where $m>i>0$. Both $L_{n}$, $L_{n+m}$ and the interpolated latent space
$L_{n+i}$ are inserted in the decoders to construct the predicted seismic data
with total CO2 injected volume $n$, $n+i$ and $n+m$, ($d_{n}^{\mathrm{pre}}$,
$d_{n+i}^{\mathrm{pre}}$ and $d_{n+m}^{\mathrm{pre}}$). The loss function of
the network is defined as

$\zeta^{g}=\sum_{j\in\\{0,i,m\\}}\sum_{n}\alpha^{\mathrm{g}}_{n+j}\left\|d_{n+j}-d_{n+j}^{\mathrm{pre}}\right\|^{2},$
(11)

where $\alpha^{\mathrm{g}}$ are the weights for seismic data at different
total CO2 injected volume for the training of interpolation network.

###  IV-C Extrapolation Network

It is worth considering how extrapolation is even possible at all for our
problems. The notion of extrapolation is the ability to make inferences that
go beyond the scope of one’s experiences. The generalization ability exhibited
by contemporary neural network algorithms is largely limited to interpolation
between data points in the training corpora so that the extrapolation can be a
much challenging task for machine learning problems [50]. However,
fundamentally different from those problems from the computer vision
community, there are unique physical laws underlying our problems. These
physical laws are, to some extent, “invariant” in the sense that the spatial
and temporal dynamics can be characterized by those laws throughout the
physical process. Thus, the extrapolation task for our problems can be enabled
by encouraging the learning of representations that reflect these “invariant”
physical laws. It has been shown in different scientific applications that
these physical laws can be effectively and implicitly captured via training
neural networks with sufficient training data [51, 52, 53].

Specific to our problems, with the interpolation network developed, we can use
it to generate more simulations to augment the original dataset for training
our extrapolation network. To better learn the complex temporal dynamics of
CO2 plume migration, we incorporate a long short-term memory (LSTM) layer
within our extrapolation network. The LSTM takes the historical latent
representations as input and predicting the future states. An illustration of
the extrapolation network is provided in Figure 3. Unlike the interpolation
network for generating the full seismic data (shown in Figure 2), the
extrapolation network focuses on the change of the data. The extrapolation
network is composed of four encoder networks and five decoder networks.
Similar to the interpolation network, both the encoders and the decoders share
their weights with other encoders and decoders. The seismic data with total
CO2 injected volume $n$, $d_{n}$ is set as the baseline and the change of the
data with total CO2 injected volume from $n$ to $n+3i$ ($\Delta{d_{(n,n)}}$,
$\Delta{d_{(n,n+i)}}$, $\Delta{d_{(n,n+2i)}}$, $\Delta{d_{(n,n+3i)}}$) are
calculated,

$\Delta{d_{(n,n+i)}}=d_{n+i}-d_{n},$ (12)

where $\Delta{d_{(n,n)}}$, $\Delta{d_{(n,n+i)}}$, $\Delta{d_{(n,n+2i)}}$,
$\Delta{d_{(n,n+3i)}}$ are the inputs for the encoder networks, and their
latent space representations $l_{n}$, $l_{n+i}$, $l_{n+2i}$ and $l_{n+3i}$ are
connected using LSTM network to predicted the latent space representation of
the change with total CO2 injected volume $n+4i$, $l_{n+4i}$. Then $l_{n+4i}$
is fit in the decoder network to generate the data change
$\Delta{d^{\mathrm{pre}}_{(n,n+4i)}}$. The loss function of our extrapolation
network is

$\zeta^{e}=\sum_{m=0}^{4}\sum_{n}\alpha^{\mathrm{e}}_{n+mi}\left\|\Delta{d_{(n,n+mi)}}-\Delta{d^{\mathrm{pre}}_{(n,n+mi)}}\right\|^{2},$
(13)

where $\alpha^{\mathrm{e}}$ is the weight for seismic data with different
total CO2 injected volume for the training of extrapolation network. In the
last step, the predicted data change $\Delta{d^{\mathrm{pre}}_{(n,n+4i)}}$ is
added to the baseline $d_{n}$ to generate the predicted data at time $n+4i$,
$d^{\mathrm{pre}}_{n+4i}$.

###  IV-D Extrapolation Network with Optical Flow Regularization

As pointed out in [40], a prominent physical phenomenon that can be observed
during the CO2 injection at Sleipner is the “velocity push-down” effect. It is
caused by the decrease of the velocity within the CO2 reservoir after the
injection of CO2 in saturating local saline water. An immediate impact of this
saturation would result in a longer traveltime for those seismic waves that
propagate through the CO2 reservoir. That in turn will result in a time-shift,
(moving downward to be more specific), in the time-domain seismic image.
During the process of CO2 injection, the velocity of the CO2 reservoir will
keep decreasing due to the increase of the saturation of CO2 within saline
water. This will result in a continuous moving-downward of the reflection
events below the CO2 reservoir, which will be consistent throughout the
injection process. That would provide a reliable constraint for extrapolation.

The technique of optical flow in Section III.C, as a method to describe the
spatial movement of an object, can both capture the migration of CO2 plume and
the downward shifting of the reflection events. Given the seismic data with
total CO2 injected volume $n$, and $n+3i$, ($d_{n}$, $d_{n+3i}$), we calculate
the observed optical flow vector $\mathbf{m_{1}}=(u_{1},v_{1},w_{1})$ from
$d_{n}$ to $d_{n+3i}$. Given the predicted data with total CO2 injected volume
$n$, $n+i$ and $n+4i$, ($d_{n}$, $d_{n+i}$, $d_{n+4i}$), we calculate the
predicted vector $\mathbf{m_{2}}=(u_{2},v_{2},w_{2})$ from
$d^{\mathrm{pre}}_{n}$ to $d^{\mathrm{pre}}_{n+4i}$ and the predicted vector
$\mathbf{m_{3}}=(u_{3},v_{3},w_{3})$ from $d^{\mathrm{pre}}_{n+i}$ to
$d^{\mathrm{pre}}_{n+4i}$. The regularization term is defined as the summation
of the included angle between $\mathbf{m_{1}}$ and $\mathbf{m_{2}}$ and the
included angle between $\mathbf{m_{1}}$ and $\mathbf{m_{3}}$ ,

$\zeta^{\mathrm{reg}}=\sum_{n}\left\|\arccos{r}_{n}+\arccos{r}_{n+1}\right\|^{2},$
(14)

where $\arccos$ is the inverse of the cosine function,

${r}_{n}=\frac{u_{1}u_{2}+v_{1}v_{2}+w_{1}w_{2}}{\sqrt{u_{1}^{2}+v_{1}^{2}+w_{1}^{2}}\sqrt{u_{2}^{2}+v_{2}^{2}+w_{2}^{2}}},$
(15)

and

${r}_{n+1}=\frac{u_{1}u_{3}+v_{1}v_{3}+w_{1}w_{3}}{\sqrt{u_{1}^{2}+v_{1}^{2}+w_{1}^{2}}\sqrt{u_{3}^{2}+v_{3}^{2}+w_{3}^{2}}}.$
(16)

Thus, the overall loss function for the extrapolation network with optical
flow regularization can be posed as

$\zeta^{\mathrm{e_{reg}}}=\zeta^{\mathrm{e}}+\alpha_{reg}\,\zeta^{\mathrm{reg}},$
(17)

where the terms of $\zeta^{\mathrm{e}}$ and $\zeta^{\mathrm{reg}}$ are
provided in Eqs. (13) and (14), respectively, and $\alpha_{reg}$ is the weight
for the regularization term. It is worth mentioning that this regularization
only constrains the direction of the seismic data change, not the amplitude of
the seismic data change.

Figure 4:  a) 2D observed data, interpolated data in the b) data domain and c)
latent space from year 1994 to 2001. For the data domain interpolation, the
interpolated data are generated by linear combination of the observed data.
For the interpolation in the latent space, the latent spaces of the observed
data are combined linearly and then fed into the decoders to generate
interpolated data. Figure 5:  The zoom views of a) 2D observed data,
interpolated data in the b) data domain and c) latent space from year 1994 to
2001.

##  V Tests

In this section, we provide numerical tests and results showing the
performance of our models in both 2D and 3D interpolation (Sections V.A & V.B)
and extrapolation (Sections V.C & V.D) scenarios. A side-by-side comparison to
two baseline models (linear model and neural network model without
regularization) is also provided.

###  V-A 2D Sleipner Interpolation

To test the 2D interpolation method, we select a 2D slice in the 3D seismic
section and collect the seismic data of this slice every monitoring year.
Among all acquired real seismic data over 8 years, a random selection of 3
subsets out of 8 existing datasets results in a total of 56 combinations
(i.e., $C_{8}^{3}$). Each subset can be used as a training sample for our
interpolation model in Figure 2. In particular, three sequential seismic data
are selected and used as the inputs (i.e., $d_{n}$, $d_{n+i}$ and $d_{n+m}$).
We use 45 of the subsets as the training set and the rest are held out for
validation. The weights in Eq. (11) ( $\alpha^{\mathrm{g}}_{n}$,
$\alpha^{\mathrm{g}}_{n+i}$ and $\alpha^{\mathrm{g}}_{n+m}$) are set as (0.2,
0.6, 0.2).

After 1,200 epochs of training, the mean-square-error (MSE) loss value
decreases to 0.06 for the training set and 0.09 for the validation set. An
example of the interpolated data is shown in Figure 4, where we provide three
real seismic datasets as the ground truth for comparison. We also provide the
synthetic seismic data generated using linear interpolation based on pixels
(namely, data domain), i.e.,

$d_{n+i}=d_{n}+\frac{I_{n+m}-I_{n+i}}{I_{n+m}-I_{n}}*(d_{n+m}-d_{n}).$ (18)

With the high fidelity synthetic seismic data being generated using our model,
we are able to “fill the gap” between each year when no real seismic data are
acquired. By looking at the zoomed view of the interpolated data in Figure 5c,
we can see the expansion of reflectors caused by the growth of CO2 plume. In
comparison, only the cross-dissolve transition can be seen with the
interpolation in the data domain in Figure 5b.

Figure 6:  a) 3D observed data and b) interpolated data in the latent space
from year 1994 to 1999. Figure 7:  The comparisons of the observed data, 2D
interpolated data, 3D interpolated data with injection volume of a) 0 Mt and
b) 4.2 Mt.

###  V-B 3D Sleipner Interpolation

We slice the 3D Sleipner data in a group of 2D seismic images along the
crossline direction, where each 2D seismic image can be regarded as a channel
image. Since there are 144 grids in the crossline direction of the 3D Sleipner
data, we have 144 channels for the input of the interpolation network. Due to
the limitation of the computational cost, 2D convolution filters are used
instead of the 3D convolution filters. The interpolation network structure is
the same as the 2D case except that the number of the channels in the input
and output layers are changed.

After 1,200 epochs, the MSE loss decreases to 0.35 for the training set and
0.45 for the validation set, both of them are higher than the losses in the 2D
case. Examples of interpolated data are shown in Figure 6 and comparisons
between the observed data, 2D interpolated data and 3D interpolated data are
shown in Figure 7, where some details are missing in the 3D interpolated data
that the red arrows point. It means the results are affected by the
computational cost problem. To achieve similar accuracy as the 2D case, 3D
convolution filter, more training epochs and a larger neural network with more
parameters are required.

Figure 8:  An example of optical flow: a) the 2D observed data in 1999, b) the
2D observed data in 2001 and c) the calculated optical flow using the data
from 1999 to 2001.  Figure 9:  Workflow to predict data between 10.0 and 12.0
Mt from data 0 between 9.9 Mt using extrapolation network  Figure 10:  a) The
observed data, b) by linear extrapolation, by extrapolation network c) without
regularization and d) with regularization from 9.6 to 12.0 Mt  Figure 11:  The
zoom view of a) observed data change, predicted data change b) by linear
extrapolation, by extrapolation network c) without regularization and d) with
regularization from 9.6 to 12.0 Mt

###  V-C 2D Sleipner Extrapolation

Only eight seismic datasets were collected in the 16 years of monitoring in
the Sleipner area which is not enough for the training of the extrapolation
network. To augment the dataset, we use the interpolation network to generate
121 seismic datasets every 0.1 Mt from 0 to 12 Mt. We assume that only the
data from 0 to 9.9 Mt are known and these data are used in the training of the
extrapolation network. The data from 10.0 to 12.0 Mt are used as the
validation sets. The change of total CO2 injected volume with each step in the
extrapolation network, $i$, is set as 0.1 Mt.

During training, we set [$\alpha^{\mathrm{e}}_{n}$,
$\alpha^{\mathrm{e}}_{n+i}$, $\alpha^{\mathrm{e}}_{n+2i}$,
$\alpha^{\mathrm{e}}_{n+3i}$, $\alpha^{\mathrm{e}}_{n+4i}$] as (0.4, 0.4, 0.4,
0.4, 1) and the optical flow weight $\alpha_{reg}$ as 0.1. 2D optical flow is
applied to calculate the shift of the seismic data caused by the injection of
CO2. An example is shown in Figure 8. The expanding of the reflection event in
the yellow box and the “push-down” effect of the events in the red boxes are
clearly detected by the optical flow.

After training for 1,200 epochs, the training loss decreases to 0.31 and the
validation loss decreases to 0.50. The workflow to predict data between 10.0
and 12.0 Mt from data 0 between 9.9 Mt is shown in Figure 9. The data from 9.6
to 9.9 Mt are fed into the extrapolation network to predict the data with 10.0
Mt. Then the predicted data with 10 Mt is combined with the data from 9.7 to
9.9 Mt to predict the data with 10.1 Mt and the rest can be done in the same
manner. As the injection volume increases, the predicted data will move
further away from the data domain we know from 0 to 9.9 Mt. As a result, the
data from the prediction will be more and more unreliable as the injection
volume increases. The predicted data from 10.0 to 12.0 Mt is shown in Figure
10 and the MSE and the mean absolute error (MAE) are shown in Table
LABEL:tab:with_reg, the MSE and MAE losses increase with the increasing of the
total CO2 injected volume.

Since the velocity change is small with a large amount of CO2 trapped in the
reservoir, the change in the seismic data is subtle when the CO2 injected
volume is over 10 Mt as shown in Figure 10. Figure 11 shows the zoom view of
the data change in the CO2 plume area from 9.6 to 12.0 Mt. The predicted data
generated by linear extrapolation in the data domain, extrapolation network
without and with optical flow regularization are compared. There are five
bright reflection events where the red arrows point in the real data. Their
amplitudes increase with the increase of total CO2 injected volume and the
changing rate varies with time. For the data generated by linear
extrapolation, the amplitudes of the events increase at a constant rate. As a
result, the data change becomes very large when the total injection volume is
12 Mt in Figure 11b, which is far away from the real data. In Figure 11c, the
predicted data change without regularization only detects three events where
the arrows point. The predicted data change is distorted when the total
injection volume is larger than 11.5 Mt. For the predicted data with
regularization in Figure 11d, all five events are predicted correctly. The
trend of the increasing amplitude and expanding of the events remains the same
even when the data is far away from the known data domain. As a result, both
the MSE and MAE losses of the predicted data with regularization are less than
the losses of the data without regularization and the data with linear
extrapolation in table LABEL:tab:with_reg.

###  V-D 3D Sleipner Extrapolation

For the 3D Sleipner data, we use the interpolation network to generate 24
seismic datasets every 0.5 Mt from 0 to 12 Mt. The data from 0 to 9.5 Mt are
used to the train the extrapolation network. The data from 10.0 to 12.0 Mt are
used as the validation set. The change of total CO2 injected volume with each
step in the extrapolation network, $i$, is 0.5 Mt. The weights
[$\alpha^{\mathrm{e}}_{n}$, $\alpha^{\mathrm{e}}_{n+i}$,
$\alpha^{\mathrm{e}}_{n+2i}$, $\alpha^{\mathrm{e}}_{n+3i}$,
$\alpha^{\mathrm{e}}_{n+4i}$] are (0.4, 0.4 ,0.4, 0.4, 1) and the optical flow
weight $\alpha_{reg}$ is 0.1.

After 600 epochs, the training loss decreases to 0.65 and the validation loss
decreases to 0.83. We use the data from 8.0 to 9.5 Mt to predict the data from
10.0 to 12.0 Mt. The losses between the 3D extrapolated data and the true data
are given in the table LABEL:tab:with_reg3D. Similar to the 2D extrapolated
data, the losses increase with the total CO2 injected volume. Moreover, they
increase much faster than the 2D cases since the change of total CO2 injected
volume in each step is much bigger. Examples of the slices in the 3D
extrapolated data from 10.0 Mt to 10.5 Mt and the data change from 8.0 Mt to
10.5 Mt are shown in Figure 12. By looking at the areas that the red arrows
points in the data change, we can see the data change generated by the
extrapolation network with regularization is closer to the true data than the
predicted data change generated by the extrapolation network without
regularization. Although the reflection events are preserved in for the
predicted data change generated linear extrapolation, cycle-skipping [54]
occurs with these events due to the “velocity push-down” effect, which causes
large MSE and MAE losses. The MSE and MAE losses of the predicted data are
given in Table LABEL:tab:with_reg3D, in which the data change predicted by the
extrapolation network with regularization has the lowest losses.

Figure 12:  2D slices of the 3D data, the data change and the data change zoom view in the yellow box of a) observed data, predicted data by b) linear extrapolation, c) without regularization and d) with regularization from 8.0 to 10.5 Mt. TABLE II: The errors in MSE and MAE between the 2D extrapolated data and observed data |  CO2 Injected Volume |  10.0 Mt |  10.5 Mt |  11.0 Mt |  11.5 Mt |  12.0 Mt
---|---|---|---|---|---|---
Predicted Data Linear Extrapolation | MSE | 0.1419 | 0.4994 | 1.1370 | 1.9074 | 2.8454
MAE | 0.1594 | 0.3191 | 0.5307 | 0.7177 | 0.8843
Predicted Data Without Regularization | MSE | 0.0529 | 0.1005 | 0.1562 | 0.1725 | 0.2178
MAE | 0.1486 | 0.2023 | 0.2751 | 0.3420 | 0.4229
Predicted Data With Regularization | MSE | 0.0488 | 0.0941 | 0.1351 | 0.1537 | 0.1621
MAE | 0.1475 | 0.1970 | 0.2569 | 0.3085 | 0.3509
TABLE III: The errors in MSE and MAE between the 3D extrapolated data and observed data |  CO2 Injected Volume |  10.0 Mt |  10.5 Mt |  11.0 Mt |  11.5 Mt |  12.0 Mt
---|---|---|---|---|---|---
Predicted Data Linear Extrapolation | MSE | 0.0475 | 0.2065 | 0.5200 | 0.9755 | 1.4381
MAE | 0.1079 | 0.2322 | 0.3751 | 0.5251 | 0.6421
Predicted Data Without Regularization | MSE | 0.1028 | 0.2064 | 0.5095 | 0.7933 | 1.1066
MAE | 0.1500 | 0.2589 | 0.4202 | 0.5418 | 0.6231
Predicted Data With Regularization | MSE | 0.0925 | 0.1338 | 0.3089 | 0.5079 | 1.0611
MAE | 0.1521 | 0.2152 | 0.3353 | 0.4401 | 0.6087

##  VI Discussions and Future Work

Expert Evaluation. Two error metrics (i.e., MSE and MAE) have been used to
evaluate the quality of our models in Tables LABEL:tab:with_reg and
LABEL:tab:with_reg3D. They are both pixel-based, which may be insufficient for
assessing structured data such as images [16]. To overcome potential biased
assessment based on MSE and MAE metrics, we provide an expert evaluation,
which will add another dimension of assessment. This is particularly important
for our work since the resulting images will eventually be presented to the
field operators for monitoring purposes. From that perspective, the standard
for the evaluation would be to pass the domain expert with our generated
samples. To rule out any subjective evaluation, we conduct a blind survey
completed by 20 domain experts to evaluate the quality of our interpolation
and extrapolation data. The participants’ number of working years in the
seismic imaging area are listed in Figure 13a, which can be used as an
indication of their level of expertise. The survey can be found in our online
Google form
111https://docs.google.com/forms/d/e/1FAIpQLSfEZVvuXHcIub_mVzmLGEIOA4ZXkFWL_Tz8sXwwG1nqRrcukA/viewform.
In this survey, we design five questions regarding specific geologic features
of overburden, CO2 plume, reflection below CO2 plume, noise, and overall
quality. We pick those features since they have been utilized to characterize
CO2 migration in various work [55, 40]. Four different seismic images
belonging to four different datasets are used in the survey including, real
data, 2D interpolated data, 3D interpolated data, and 2D extrapolated data. We
then ask the experts to grade (on a scale of 4 points) the quality of the
images based on the aforementioned geologic features without telling the
participants in advance whether the data is generated or not. Average scores
of the expert evaluation are shown in Figure 13b.

For the overburden area, the score of true data is 4, which is much higher
than the scores of 2D interpolated data, 3D interpolated data and 2D
extrapolated data. This indicates the overburden area of the generated data
needs further improvements. However, it is interesting to notice that for the
rest of the geologic features, the scores are all very close. For example, the
scores for the overall quality of the real data, 2D interpolated data, 3D
interpolated data and 2D extrapolated data are [3.05, 3.15, 2.95, 3.00], which
means the quality of the generated data is close enough to the real data. Even
the experts cannot distinguish them. We believe that shows the excellent
performance of our models.

Figure 13:  a) Working years of the experts in the seismic imaging area. b)
Average scores of the expert evaluation on the quality of the data.

Space and Time Trade-off with a More Efficient 3D Implementation. In our
current implementation of the 3D model, our focus is to demonstrate the
effectiveness of our idea in both interpolation and extrapolation for 2D/3D
scenarios. The performance for 3D scenarios is indeed worse than 2D because of
the limitations of the computational resources in both memory and efficiency.
When extending the model from 2D model to 3D, the major technical barriers
would be the high computational training cost and the large memory
consumption. Although the additional computational improvement for 3D
scenarios is not the focus of this manuscript, we have already identified
numerical solutions to overcome the limitations of computational resources.
Particularly, we will leverage group convolution techniques to reduce the
global convolution operations on the channel dimension, which has been proved
to be effective in reducing the overall cost [56]. To reduce the memory
consumption, we will trade storage for computation. It is well known that most
of the memory consumption of a neural network models result from a large
number of intermediate activations. Instead of saving all the activations, we
will only store some activations and recompute those in-between when needed.
This is described as in the “Invertible Network” [57].

Known Knowns versus Unknown Unknowns. Extrapolation, in general, is
challenging for any machine learning model including ours. However, what makes
our problem unique and different comparing to many other problems is the
existence of the underlying physics knowledge, which yields consistent spatio-
temporal dynamics such as “velocity push-down”. Our neural network models have
been proved to be effective in capturing those characteristics. That in turn
will help in the prediction of the near-future unknown data due to the
invariant nature of the underlying physics. This is what we refer to as “known
knowns”. There will also be “unknown unknowns”, which refers to uncertainty in
the data. Provided with the limited availability of the existing data, the
limitation of our model, and some additional unpredictable variances of the
problem itself, the predictability of the machine learning models would become
unattainable. That has also been reflected in our numerical testing results
(i.e., the increase of the loss values furthering out in the future). To
improve the prediction performance, one immediate solution would be to fine-
tune the model with newly acquired data.

Generative Models and Others. Deep generative models (DGMs), such as
generative adversarial network (GAN) [58, 32], variational autoencoder (VAE)
[49] and adversarial autoencoder (AAE) [59], learn lower-dimensional data
distribution explicitly or implicitly. Once the data distribution is
successfully captured, DGMs can be used to generate synthetic manipulated data
for interpolation. However, training DGMs requires a great amount of data,
which would make them inappropriate for the interpolation and extrapolation of
Sleipner data. For cases when a lot more 4D seismic datasets with much more
time samples become available, DGMs would be interesting approaches to study.
We plan to study their capability in spatio-temporal in-between image
generation in the future. In our method, we utilize optical flow to track the
temporal change of the seismic data over time and incorporate it to regularize
the loss function. We will compare the performance using other motion tracking
techniques, such as template matching [60], dynamic warping [61, 62] and
kernelized correlation filter [63]. In addition, it will be interesting to
incorporate non-seismic data in the Sleipner area, such as gravity data and
well-logging, into our proposed network to further regularize our model and
potentially improve the results. The selection of the hyper-parameters (such
as $\alpha^{\mathrm{g}}_{n}$, $\alpha^{\mathrm{g}}_{n+i}$,
$\alpha^{\mathrm{g}}_{n+m}$ used in Eq. (11) and $\alpha^{\mathrm{e}}_{n}$,
$\alpha^{\mathrm{e}}_{n+i}$, $\alpha^{\mathrm{e}}_{n+2i}$,
$\alpha^{\mathrm{e}}_{n+3i}$, $\alpha^{\mathrm{e}}_{n+4i}$ used in Eq. (13))
play an important role in producing accurate estimation. In this manuscript,
we set the parameter based on their importance to the resulting prediction.
Considering the relatively small size of our dataset, our empirical approach
to select those hyper-parameters has been demonstrated to be rather effective
via the numerical results. For a large dataset with more varying spatio-
temporal dynamics, a more heuristic optimization technique would be needed to
derive the optimal parameters. We would consider that as a future direction.

##  VII Conclusions

In this paper, we develop two neural-network-based methods for spatio-temporal
in-between image interpolation and extrapolation. To capture the spatial
representation, our models leverage the deep autoencoder to obtain
semantically meaningful latent variables, which can be combined linearly to
generate new features for interpolation/extrapolation tasks. Temporal dynamics
are also important in designing our models. To enhance the temporal
correlation of data, we utilize both the LSTM structure and optical flow
regularization in our neural network models. We validate the performance of
our models using 2D/3D post-stack seismic imaging data acquired at the
Sleipner CO2 sequestration field. Through both numerical comparisons to
baseline methods and experts evaluation, we show our models yield high
fidelity seismic imaging data that would allow in-situ reservoir monitoring
and forecasting, and therefore yielding the great potential for CO2
sequestration and storage.

##  VIII Acknowledgements

This work was co-funded by the U.S. DOE Office of Fossil Energy’s Carbon
Storage program and by the Laboratory Directed Research and Development
program of LANL under project numbers 20210542MFR and 20200061DR.

## References

  * [1] J. Ma, L. Lin, H. Wng, M. Tan, S. Cui, Y. Zhang, Z. Qu, L. Jia, and S. Zhang, “Geophysical monitoring technology for CO2 sequestration,” _Applied Geophysics_ , vol. 13, no. 2, pp. 288–306, 2016.
  * [2] K. Koster, P. Gabriels, M. Hartung, J. Verbeek, G. Deinum, and R. Staples, “Time-lapse seismic surveys in the North Sea and their business impact,” _The Leading Edge_ , vol. 19, no. 3, pp. 286–293, 2000.
  * [3] R. Chadwick, P. Zweigel, U. Gregersen, G. Kirby, S. Holloway, and P. Johannessen, “Geological reservoir characterization of a CO2 storage site: The Utsira sand, Sleipner, northern North Sea,” _Energy_ , vol. 29, no. 9-10, pp. 1371–1381, 2004.
  * [4] A. Chadwick, G. Williams, N. Delepine, V. Clochard, K. Labat, S. Sturton, M.-L. Buddensiek, M. Dillen, M. Nickel, A. L. Lima _et al._ , “Quantitative analysis of time-lapse seismic monitoring data at the Sleipner CO2 storage operation,” _The Leading Edge_ , vol. 29, no. 2, pp. 170–177, 2010\.
  * [5] P. Cristovao, H. Nakada, Y. Tanimura, and H. Asoh, “Generating in-between images through learned latent space representation using variational autoencoders,” _IEEE Access_ , vol. 8, pp. 149 456–149 467, 2020.
  * [6] H. Jiang, D. Sun, V. Jampani, M. Yang, E. Miller, and J. Kautz, “Super SloMo: High quality estimation of multiple intermediate frames for video interpolation,” in _IEEE/CVF Conf. Comput. Vis. Pattern Recognit._ , 2018, pp. 9000–9008.
  * [7] L. P. Zuckerman, E. Naor, G. Pisha, S. Bagon, and M. Irani, “Across scales and across dimensions: Temporal super-resolution using deep internal learning,” in _European Conference on Computer Vision_. Springer, 2020, pp. 52–68.
  * [8] A. Oring, Z. Yakhini, and Y. Hel-Or, “Faithful autoencoder interpolation by shaping the latent space,” _arXiv preprint arXiv:2008.01487_ , 2020.
  * [9] D. W. Vasco, A. Datta-Gupta, R. Behrens, P. Condon, and J. Rickett, “Seismic imaging of reservoir flow properties: Time-lapse amplitude changes,” _Geophysics_ , vol. 69, no. 6, pp. 1425–1442, 2004.
  * [10] C. MacBeth, H. Amini, and S. Izadian, “Methods of measurement for 4D seismic post-stack time shifts,” _Geophysical Prospecting_ , vol. 68, no. 9, pp. 2637–2664, 2020.
  * [11] D. E. Lumley, “Business and technology challenges for 4D seismic reservoir monitoring,” _The Leading Edge_ , vol. 23, no. 11, pp. 1166–1168, 2004.
  * [12] A.-K. Furre, O. Eiken, H. Alnes, J. N. Vevatne, and A. F. Kiær, “20 years of monitoring CO2-injection at Sleipner,” _Energy procedia_ , vol. 114, pp. 3916–3926, 2017.
  * [13] R. Chadwick, R. Arts, and O. Eiken, “4D seismic quantification of a growing CO2 plume at Sleipner, North Sea,” in _Geological Society, London, Petroleum Geology Conference series_ , vol. 6, no. 1. Geological Society of London, 2005, pp. 1385–1399.
  * [14] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [15] Z. Li and A. Abubakar, “Complete sequence stratigraphy from seismic optical flow without human labeling,” in _SEG Technical Program Expanded Abstracts_. Society of Exploration Geophysicists, 2020, pp. 1248–1252.
  * [16] R. Zhang, P. Isola, A. A. Efros, E. Shechtman, and O. Wang, “The unreasonable effectiveness of deep features as a perceptual metric,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 586–595.
  * [17] S. Spitz, “Seismic trace interpolation in the FX domain,” _Geophysics_ , vol. 56, no. 6, pp. 785–794, 1991.
  * [18] Y. Jia and J. Ma, “What can machine learning do for seismic data processing? An interpolation application,” _Geophysics_ , vol. 82, no. 3, pp. V163–V177, 2017.
  * [19] K. Lu and S. Feng, “Auto-windowed super-virtual interferometry via machine learning: A strategy of first-arrival traveltime automatic picking for noisy seismic data,” in _SEG 2018 Workshop: SEG Maximizing Asset Value Through Artificial Intelligence and Machine Learning, Beijing, China, 17-19 September 2018_. Society of Exploration Geophysicists and the Chinese Geophysical Society, 2018, pp. 10–14.
  * [20] S. Mandelli, V. Lipari, P. Bestagini, and S. Tubaro, “Interpolation and denoising of seismic data using convolutional neural networks,” _arXiv preprint arXiv:1901.07927_ , 2019.
  * [21] H. Kaur, N. Pham, and S. Fomel, “Seismic data interpolation using CycleGAN,” in _SEG Technical Program Expanded Abstracts_. Society of Exploration Geophysicists, 2019, pp. 2202–2206.
  * [22] T. Liu, H. Chen, G. Hetz, and A. Datta-Gupta, “Integration of time-lapse seismic data using the onset time approach: The impact of seismic survey frequency,” _Journal of Petroleum Science and Engineering_ , vol. 189, p. 106989, 2020.
  * [23] A. F. Kiær, O. Eiken, and M. Landrø, “Calendar time interpolation of amplitude maps from 4D seismic data,” _Geophysical Prospecting_ , vol. 64, no. 2, pp. 421–430, 2016.
  * [24] E. Lindeberg, P. Zweigel, P. Bergmo, A. Ghaderi, and A. Lothe, “Prediction of CO2 distribution pattern improved by geology and reservoir simulation and verified by time lapse seismic,” _Greenhouse gas control technologies_ , vol. 372, p. 377, 2001.
  * [25] G. M. Hoversten, R. Gritto, J. Washbourne, and T. Daley, “Pressure and fluid saturation prediction in a multicomponent reservoir using combined seismic and electromagnetic imaging,” _Geophysics_ , vol. 68, no. 5, pp. 1580–1591, 2003.
  * [26] P. Upchurch, J. Gardner, G. Pleiss, R. Pless, N. Snavely, K. Bala, and K. Weinberger, “Deep feature interpolation for image content changes,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 7064–7073.
  * [27] D. Wang, A. Vincent, P. Blanchfield, , and R. Klepko, “Motion compensated frame rate up-conversion - Part II: New algorithms for frame interpolation,” _IEEE Trans. Broadcast._ , vol. 56, no. 2, pp. 142–149, 2010\.
  * [28] J. Zhai, K. Yu, J. Li, and S. Li, “A low complexity motion compensated frame interpolation method,” in _IEEE Int. Symp. Circuits Syst._ , 2005, pp. 4927–4930.
  * [29] S. Meyer, O. Wang, H. Zimmer, M. Grosse, and A. Sorkine-Hornung, “Phase-based frame interpolation for video,” in _IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR)_ , 2015, pp. 1410–1418.
  * [30] N. Wadhwa, M. Rubinstein, F. Durand, and W. Freeman, “Phase-based video motion processing,” _ACM Trans. Graph._ , vol. 32, no. 4, pp. 1–10, 2013.
  * [31] S. Meyer, A. Djelouah, B. McWilliams, A. Sorkine-Hornung, M. Gross, and C. Schroers, “PhaseNet for video frame interpolation,” in _IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR)_ , 2018, pp. 498–507.
  * [32] D. Berthelot, C. Raffel, A. Roy, and I. Goodfellow, “Understanding and improving interpolation in autoencoders via an adversarial regularizer,” _arXiv preprint arXiv:1807.07543_ , 2018.
  * [33] T. Sainburg, M. Thielk, B. Theilman, B. Migliori, and T. Gentner, “Generative adversarial interpolative autoencoding: adversarial training on latent space interpolations encourage convex latent distributions,” _arXiv preprint arXiv:1807.06650_ , 2018.
  * [34] E. Ilg, N. Mayer, T. Saikia, M. Keuper, A. Dosovitskiy, and T. Brox, “FlowNet 2.0: Evolution of optical flow estimation with deep networks,” in _IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR)_ , 2017, pp. 2462–2470.
  * [35] T. White, “Sampling generative networks,” _arXiv preprint arXiv:1609.04468_ , 2016.
  * [36] D. Lumley, “4D seismic monitoring of CO2 sequestration,” _The Leading Edge_ , vol. 29, no. 2, pp. 150–155, 2010.
  * [37] C. Sambo, C. C. Iferobi, A. A. Babasafari, S. Rezaei, and O. A. Akanni, “The role of 4D time lapse seismic technology as reservoir monitoring and surveillance tool: A comprehensive review,” _Journal of Natural Gas Science and Engineering_ , p. 103312, 2020.
  * [38] R. Behrens, P. Condon, W. Haworth, M. Bergeron, and C. Ecker, “4D seismic monitoring of water influx at bay marchand: The practical use of 4D in an imperfect world,” _SPE Reservoir Evaluation & Engineering_, vol. 5, no. 05, pp. 410–420, 2002.
  * [39] A. Ivanova, A. Kashubin, N. Juhojuntti, J. Kummerow, J. Henninges, C. Juhlin, S. Lüth, and M. Ivandic, “Monitoring and volumetric estimation of injected CO2 using 4D seismic, petrophysical data, core measurements and well logging: a case study at Ketzin, Germany,” _Geophysical Prospecting_ , vol. 60, no. 5, pp. 957–973, 2012.
  * [40] F. Boait, N. White, M. Bickle, R. Chadwick, J. Neufeld, and H. Huppert, “Spatial and temporal evolution of injected CO2 at the Sleipner Field, North sea,” _Journal of geophysical research: solid earth_ , vol. 117, no. B3, 2012.
  * [41] P. Zweigel, R. Arts, A. E. Lothe, and E. B. Lindeberg, “Reservoir geology of the Utsira formation at the first industrial-scale underground CO2 storage site (Sleipner area, North Sea),” _Geological Society, London, Special Publications_ , vol. 233, no. 1, pp. 165–180, 2004.
  * [42] A. J. Cavanagh and R. S. Haszeldine, “The Sleipner storage site: Capillary flow modeling of a layered CO2 plume requires fractured shale barriers within the Utsira formation,” _International Journal of Greenhouse Gas Control_ , vol. 21, pp. 101–112, 2014.
  * [43] O. Eiken, “Twenty years of monitoring CO2 injection,” _Geophysics and geosequestration_ , p. 209, 2019.
  * [44] S. Feng, O. Yilmaz, Y. Chen, and G. T. Schuster, “Zero-offset sections with a deblurring filter in the time domain,” _Geophysics_ , vol. 84, no. 4, pp. S239–S249, 2019.
  * [45] R. Chadwick, G. Williams, and I. Falcon-Suarez, “Forensic mapping of seismic velocity heterogeneity in a CO2 layer at the Sleipner CO2 storage operation, North sea, using time-lapse seismics,” _International Journal of Greenhouse Gas Control_ , vol. 90, p. 102793, 2019.
  * [46] “Utsira TWT map,” https://co2datashare.org/landing_pages/sleipner-4d-seismic-dataset/static/SleipnerSeismicMap.png, 2020\.
  * [47] D. Bank, N. Koenigstein, and R. Giryes, “Autoencoders,” _arXiv preprint arXiv:2003.05991_ , 2020.
  * [48] Q. Zhang, “RTM angle gathers and specular filter (SF) RTM using optical flow,” in _SEG Technical Program Expanded Abstracts_. Society of Exploration Geophysicists, 2014, pp. 3816–3820.
  * [49] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” _stat_ , vol. 1050, p. 10, 2014.
  * [50] C. Bishop, _Pattern Recognition and Machine Learning_. Springer Science & Business Media, 2006.
  * [51] L. Sun, H. Gao, S. Pan, and J.-X. Wang, “Surrogate modeling for fluid flows based on physics-constrained deep learning without simulation data,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 361, p. 112732, 2020.
  * [52] R. Rojas-Gomez, J. Yang, Y. Lin, J. Theiler, and B. Wohlberg, “Physics-consistent data-driven waveform inversion with adaptive data augmentation,” _IEEE Geoscience and Remote Sensing Letters_ , 2020.
  * [53] S. Feng, L. Passone, and G. T. Schuster, “Superpixel-based convolutional neural network for georeferencing the drone images,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 14, pp. 3361–3372, 2021.
  * [54] Y. Chen, Z. Feng, L. Fu, A. AlTheyab, S. Feng, and G. Schuster, “Multiscale reflection phase inversion with migration deconvolution,” _Geophysics_ , vol. 85, no. 1, pp. R55–R73, 2020.
  * [55] R. Chadwick, B. Marchant, and G. Williams, “CO2 storage monitoring: leakage detection and measurement in subsurface volumes from 3D seismic data at Sleipner,” _Energy Procedia_ , vol. 63, pp. 4224–4239, 2014.
  * [56] A. G. Howard, M. Zhu, B. Chen, D. Kalenichenko, W. Wang, T. Weyand, M. Andreetto, and H. Adam, “Mobilenets: Efficient convolutional neural networks for mobile vision applications,” _arXiv preprint arXiv:1704.04861_ , 2017.
  * [57] L. Dinh, D. Krueger, and Y. Bengio, “Nice: Non-linear independent components estimation,” _arXiv preprint arXiv:1410.8516_ , 2014.
  * [58] A. Radford, L. Metz, and S. Chintala, “Unsupervised representation learning with deep convolutional generative adversarial networks,” _arXiv preprint arXiv:1511.06434_ , 2015.
  * [59] A. Makhzani, J. Shlens, N. Jaitly, I. Goodfellow, and B. Frey, “Adversarial autoencoders,” _arXiv preprint arXiv:1511.05644_ , 2015.
  * [60] R. Brunelli, _Template matching techniques in computer vision: theory and practice_. John Wiley & Sons, 2009.
  * [61] D. Hale, “Dynamic warping of seismic images,” _Geophysics_ , vol. 78, no. 2, pp. S105–S115, 2013.
  * [62] S. Feng and G. T. Schuster, “Transmission+ reflection anisotropic wave-equation traveltime and waveform inversion,” _Geophysical Prospecting_ , vol. 67, no. 2, pp. 423–442, 2019.
  * [63] J. F. Henriques, R. Caseiro, P. Martins, and J. Batista, “High-speed tracking with kernelized correlation filters,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 37, no. 3, pp. 583–596, 2014.


# Quantum transport of strongly interacting fermions in one dimension at far-
out-of-equilibrium

Jie Zou  State Key Laboratory of Surface Physics, Institute of Nanoelectronics
and Quantum Computing, and Department of Physics, Fudan University, Shanghai
200433, China  Xiaopeng Li<EMAIL_ADDRESS>State Key Laboratory of
Surface Physics, Institute of Nanoelectronics and Quantum Computing, and
Department of Physics, Fudan University, Shanghai 200433, China  Shanghai Qi
Zhi Institute, AI Tower, Xuhui District, Shanghai 200232, China

###### Abstract

In the study of quantum transport, much has been known for dynamics near
thermal equilibrium. However, quantum transport far away from equilibrium is
much less well understood—the linear response approximation does not hold for
physics far-out-of-equilibrium in general. In this work, motivated by recent
cold atom experiments on probing quantum many-body dynamics of a one-
dimensional XXZ spin chain, we study the strong interaction limit of the one-
dimensional spinless fermion model, which is dual to the XXZ spin chain. We
develop a highly efficient computation algorithm for simulating the non-
equilibrium dynamics of this system exactly, and examine the non-equilibrium
dynamics starting from a density modulation quantum state. We find ballistic
transport in this strongly correlated setting, and show a plane-wave
description emerges at long-time evolution. We also observe sharp distinction
between transport velocities in short and long times as induced by interaction
effects, and provide a quantitative interpretation for the long-time transport
velocity.

Introduction.— Transport in low dimensional quantum systems [1] has been
attracting continuous research interests from both theoretical and
experimental perspectives, due to both the existence of solvable models [2]
and intrinsic strong correlation effects [1, 3]. Although wavefunctions are
exactly solvable via Bethe ansatz (BA) [4] in certain one-dimensional models,
it is challenging to compute dynamical observables such as transport
properties with analytic methods. Novel phenomenon unveiled by recent studies
still demand advanced numerical techniques to explain [5, 6, 7].

The spin-$1/2$ XXZ chain in one dimension is a representative spin model, with
the Hamiltonian

$H=J\sum_{i=1}^{L}(S_{i}^{x}S_{i+1}^{x}+S_{i}^{y}S_{i+1}^{y}+\Delta
S_{i}^{z}S_{i+1}^{z}),$ (1)

where $L$ is the system size, $J$ denotes the hopping strength and $\Delta$
labels the spin anisotropy. Despite its simplicity, this model exhibits
several interesting transport properties worth theoretical investigation. New
concepts in nonequilibrium dynamics in closed quantum systems [8, 9],
including generalized Gibbs ensemble [10], can also be studied in this model
with a quantum quench setup [11]. Another reason is from the perspective of
phenomenological transport laws. The existence of anomalous transport at
finite temperatures, which does not obey Fourier’s law [12], has been proved
via linear-response theory [13] in this model.

Previous studies on one-dimensional quantum spin transport in the linear
response regime have provided important insights about the correlated quantum
dynamics. It has been proved that finite-temperature transport is ballistic
for $|\Delta|<1$ (gapless phase) at zero magnetization
($m_{z}=\sum_{i}\braket{s_{i}^{z}}/L=0$) [14] or for generic $\Delta$ at
finite magnetization [15, 13]. For the former, recent studies suggest [16, 17]
a diffusive transport in gapped phase ($|\Delta|>1$) and predict [18] a
superdiffusive transport in the Heisenberg limit ($|\Delta|\rightarrow 1$)
when $m_{z}=0$. For the latter, Mazur inequality [19] gives a finite lower
bound [15, 13] for the spin Drude weight $D_{W}^{(s)}$ [20], a criterion of
ballistic transport. The theoretical study at high temperature [13] finds

$D_{W}^{(s)}\geq\frac{J^{2}\Delta^{2}m_{z}^{2}}{4T}\frac{1-m_{z}^{2}}{1+\Delta^{2}(2+2m_{z}^{2})}>0.$
(2)

In contrast to the linear-response regime, dynamics far-out-of-equilibrium is
less well studied. The results obtained via linear-response theory may not
hold at far-out-of-equilibrium. Numerical methods, such as the commonly used
density matrix renormalization group (DMRG) algorithm [21], tend to be limited
in simulating long time quantum dynamics due to the rapid entanglement growth.
Finding solvable models or efficient algorithms is thus important for
investigating far-out-of-equilibrium dynamics.

The study of far-out-of-equilibrium quantum transport becomes more desirable
with recent developments in cold atom experiments. Experiments based on such
clean and controllable closed quantum systems [22] are naturally suitable to
study the nonequilibrium dynamics, stimulating much recent research efforts in
this direction [23, 5, 7]. Notably, a recent experiment [7] demonstrates a
subdiffusive transport for spin helix state in a gapped phase at zero
magnetization for the anisotropic Heisenberg model, instead of a diffusive
behavior as expected from linear-response theory analysis [16]. The experiment
on the isotropic case also shows a violation of the linear-response
approximation when considering high-energy-density spin spiral states [5].
Modeling and understanding the far-out-of-equilibrium quantum transport of the
one dimensional spin chain demands further theoretical investigation.

Figure 1: (a) Workflow of the algorithm solving the strongly interacting
model. The dynamics can be efficiently computed in logical basis by mapping
the interacting model to a free one. The physical observables are computed by
mapping back to interacting basis [24, 25], for which the key is to sample the
Slater determinant wavefunction efficiently. Here we implement the recently
developed fast fermion sampling algorithm [26], in order to simulate large
system sizes. (b) Illustration of the density modulation state. In this state,
the particles are distributed in a periodic way with a period, $\lambda$, and
are separated by a distance, $R$.

In this work, we focus on a one-dimensional spinless fermion model, which is
equivalent to the spin-$1/2$ XXZ model. We take a strong interaction limit, in
which physical observables can be efficiently calculated with our fast fermion
sampling algorithm [24, 25, 26] However, due to the infinitely large
repulsion, we can only deal with cases away from half-filling (zero
magnetization). We choose the density modulation state (DMS), a highly excited
state, to check whether linear response results still hold in the far-out-of-
equilibrium setting. The setup of DMS state is demonstrated in Fig. 1(b).
Since the average energy of DMS corresponds to that of an infinite-temperature
ensemble, the lower bound of spin Drude weight (Eq. (2)) vanishes. In spite of
this vanishing lower bound, our results show that the system still exhibits
ballistic transport. Furthermore, we provide a plane-wave picture in
describing the long-time dynamics of this highly excited state in the strongly
interacting model. At last, we observe different transport velocities in short
and long times caused by interaction effects. The transport velocity in the
long-time limit is described by our plane-wave picture.

Model and Methods.— We consider a strongly interacting fermionic model on a
one-dimensional lattice, with particle number $N$ and the site index
$i\in[1,L]$. The quantum many-body dynamics is described by a Hamiltonian,

$H=J\sum_{i=1}^{L_{1}}(c^{\dagger}_{i}c_{i+1}+h.c.)+U\sum_{i=1}^{L_{1}}\sum_{k=1}^{R}n_{i}n_{i+k},$
(3)

with $i$ the lattice site index, $c_{i}$ the fermionic annihilation operator,
and $n_{i}=c_{i}^{\dagger}c_{i}$ the occupation number. Here we take the
strongly interacting limit $U\rightarrow+\infty$ and $J$ is the tunneling
strength. The interaction range $R$ and the filling $\rho=N/L$ are both
tunable parameters in our model. In this work, we choose an open boundary
condition and set the lattice constant and tunneling to be the length and
energy units, respectively. In fact, this fermionic model is equivalent to a
spin-1/2 XXZ model through a Jordan-Wigner transformation [4]. Notably, all
the dynamics are suppressed in this strong interaction limit when the system
is at half-filling, or equivalently at zero magnetization in the corresponding
spin model. Therefore, in this work, we carry out our theoretical study away
from half-filling.

To overcome the exponential complexity of the interacting many-body system,
here we adapt an efficient sampling algorithm [25, 26] to compute the time
evolution, with an $O(M_{s}LN^{2})$ complexity and an error $\delta$
proportional to the inverse square root of $M_{s}$, the number of random
samples. The illustration of this method is shown in Fig. 1(a). Since we take
the infinite interaction limit, those states with two particles in a
neighborhood of size $R$ are prohibited. Taking this constraint into account,
we can map our model to an exactly solvable free fermion case [24, 25]. This
is realized by considering one occupied site and $R$ empty sites on its left
as a a composite fermion. Through this mapping, the particle number remains
the same as the original model and the system size is reduced down to
$L_{2}=L_{1}-R(N-1)$. Hereafter, we refer to the model before (after) the
mapping as the one in physical (logical) basis. Furthermore, an advanced
fermion sampling method benefits us with less time cost when solving the free
fermion problem. Contrary to the commonly used density matrix renormalization
group (DMRG) [21] method, which only allows an investigation on systems up to
hundreds of sites in short and intermediate times, our algorithm can solve the
dynamics exactly regardless of the evolution time, for a large system with up
to thousands of sites.

Taking the interacting fermion model, we consider the time evolution of a
density modulation state, which are analogous to spin helix states
investigated in recent experiments [5, 7] through Jordan Wigner
transformation. It is a highly excited state since the expectation value of
its energy is equal to zero, corresponding to the infinite temperature
ensemble average. The preparation of this state is shown in Fig. 1(b).
Fermions are arranged in a periodic structure with $\lambda$ the wavelength,
and thus the particle number in a period is determined for a given filling
$\rho$. Within a period, these particles are aligned in a most compact
pattern, that is, the distance of two adjacent particles is set to be the
interaction range $R$. At last, a suitable system size is also necessary to
minimize the finite size effect and the computation cost, which we set
$L=4\lambda$ in this work.

Figure 2: Relaxation of the density modulation state. The interaction range
here is fixed to $R=1$. (a, b), the normalized contrast, $c(t)$ in Eq. (5),
with different wavevectors ($Q$), at filling $\rho=1/3$. The evolution time
$t$ is rescaled by the wavelength $\lambda$ and its square $\lambda^{2}$ in
(a) and (b), respectively. (c, d), the time evolution of the normalized
contrast at a different filling, $\rho=7/15$. (e) Long-time oscillation
frequency ($\omega$) of the contrast as a function of the wavevector ($Q$).
The numerical data points are well fitted by a linear function
$\omega=v_{Q}Q$, with $v_{Q}$ the dispersion velocity. The slopes $v_{Q}$ are
2.84, 3.21, 3.49, 3.69, and 3.71 for filling factors $\rho=1/3$, $2/5$, $4/9$,
$9/15$, and $10/21$. (f), the dependence of $v_{Q}$ on the filling. A
nontrivial linear relation between the inverse of dispersion velocity and the
filling factor is uncovered. We find $1/v_{Q}\approx 0.52\times(1-\rho)$ by an
empirical fitting.

Relaxation.— The relaxation dynamics is characterized by examination of the
density modulation amplitude. Due to the periodicity of density modulation
states, the contrast of density modulation is defined through a Fourier
transform

$C(t)=\frac{2}{L}\sum_{i=1}^{L}\braket{n_{i}(t)}\cos(Qi),$ (4)

where $Q=2\pi/\lambda$ is the wavevector. For comparison, we focus on the
normalized contrast

$c(t)=C(t)/C(0).$ (5)

This physical observable is also commonly used to characterize a spin helix
state in recent experiments [7].

Figure 3:  The dependence of density modulation relaxation dynamics on the
interaction range. (a, b, c), the normalized contrast $c(t)$ (Eq. (5)) with
interaction ranges $R=1$, $R=2$, and $R=3$. We investigate the relaxation
dynamics for a broad range of wavevector ($Q$) from $0.09$ to $0.16$. The
evolution time $t$ is rescaled by the wavelength $\lambda$. (d), long-time
oscillation frequency of the contrast. The numerical data points are all
captured by the linear relation $\omega=v_{Q}Q$. The dispersion velocities
$v_{Q}$, are 2.44, 3.28, and 4.77, for the interaction range $R=1$, $2$, and
$3$, respectively. Its inset shows the linear relation between the inverse of
the dispersion velocity $v_{Q}$ and the interaction range $R$. The filling
$\rho$ is fixed to $1/5$ here.

In the linear response regime, whether a system is half-filled plays a key
role in its transport property [13, 14]. One of the main purposes of our work
is to investigate the strong interaction effects on the quantum transport. We
set the filling $\rho=\frac{m}{2m+1}\,(m=1,2,...)$ so that we can gradually
probe the half-filling limit. The results are shown in Fig. 2. We observe a
perfect data collapse in the density wave contrast ($C(t)$ in Eq. (4)), with
the evolution time ($t$) rescaled by the wavelength $\lambda$. This data
collapse holds even for a filling factor quite close to $1/2$ ($\rho=7/15$).
This data-collapse behavior is absent if the time $t$ is rescaled by the
square of the wavelength $\lambda^{2}$ instead. These numerical results imply
universal ballistic transport in this model, despite the strong interaction
induced nonlinearity.

Since the quantum transport in this model is ballistic, we make an analogue
with plane wave propagating in free space with a certain velocity.
Surprisingly, the ballistic transport that emerges from the far-out-of-
equilibrium DMS state in our strongly interacting model is captured by a
physical picture of propagating plane-wave in free-space. For a massless plane
wave mode spreading in free space with wavevector $k$ and energy
$E=\hbar\omega$, the dispersion relation is $\omega=vk$, where $v$ is a
dispersion velocity that determines the ballistic transport velocity of the
plane wave in free space. In our model, we find the linear dispersion emerges
$\omega_{Q}=v_{Q}Q$ in the long-time oscillation dynamics of the DMS state,
for a broad range of wavevector, $Q$ (Fig. 2). Besides, the dispersion
velocity extracted from the long-time oscillation quantitatively agrees with
the spreading velocity (see Fig. 4(c)). These numerical results imply the
quantum dynamics of the strongly interacting fermions in one dimension
starting from the far-out-of-equilibrium DMS state is well captured by a
physical picture of propagating dynamics of a massless plane wave mode. The
strong interaction effects on the dynamics are reflected by the nontrivial
dependence of the dispersion velocity $v_{Q}$ on the filling $\rho$ (Fig.
2(f)), to be further elaborated in the later analysis of velocity
renormalization.

To confirm the generality of our finding, we further investigate the
relaxation dynamics at different interaction ranges ($R$). The results are
shown in Fig. 3. We observe that for $R=1$, $2$, and $3$, the relaxation
dynamics starting from DMS states with different wavelengths ($\lambda$) all
collapse on a curve when time $t$ is rescaled by $\lambda$. This confirms the
ballistic transport for the different choices of interaction ranges. We find
that the dispersion velocity extracted from late time oscillation dynamics
becomes larger as we increase the interaction range $R$ (Fig. 3(d)). This can
be attributed to that longer interaction range in our model tends to
accelerate the dynamical transport.

Figure 4: Sound velocity renormalization. (a) and (b) show the time evolution
of the density profile $\braket{n_{i}(t)}$ for free and interacting fermions,
respectively. The sound velocities are determined according to the propagation
of the wave front as marked by the green line in (a), and the blue and red
lines in (b). The free fermion dynamics is characterized by one velocity
scale, $v_{0}$, whereas the two different velocities, $v_{1}$ and $v_{2}$,
appear in the early and late stages of the interacting dynamics. In (b), we
choose a system size $L=96$, particle number $N=32$, wavelength $\lambda=24$.
The initial state is a density modulation state as shown in Fig. 1. We use the
same setup for non-interacting fermions in (a) for a fair comparison. (c),
convergence of the propagating velocities with increasing the system size. The
velocities for the non-interacting model (black triangles and stars), with
different fillings $\rho$ and particle separations within a period $D$,
converge to $v=1.96$ (black dotted line). The velocity $v_{2}$ for interacting
fermions (blue dots and red triangles) converges to the dispersion velocity
$v_{Q}$ (dashed lines), as extracted from long time oscillation of the
dynamics. In (a), (b), and (c), we choose $R=1$. (d) The ratio between $v_{0}$
and $v_{2}$ as a function of the filling $\rho$. The numerical results are
well fitted by a linear function.

Velocity Renormalization.— Velocity is a characteristic property of ballistic
transport, as well as an observable which can be measured directly in
experiments. In our numerical simulation here, we determine the sound velocity
by extracting the slope of the stripe pattern in the density profile
$\braket{n_{i}(t)}$. A typical time evolution of $\braket{n_{i}(t)}$ is shown
in Fig. 4(a) ((b)), considering the density modulation state in the free
(strongly interacting) model. The free fermion velocity $v_{0}$ and the short-
time (long-time) velocity $v_{1}$($v_{2}$) in the interacting case is also
determined by fitting the stripe pattern in Fig. 4(a) and (b). The strongly
interacting model shows different behaviors at early and late time evolution.
This two-stage dynamical feature as originated from strong interaction effects
makes the system distinctive from noninteracting fermions. Moreover, even at
late time where the ballistic transport of the strongly interacting model
looks similar to the non-interacting case, the propagating velocity is
actually strongly renormalized by the interaction.

To gain more insight into the renormalization process, we quantitatively study
the velocities mentioned above. As system size grows, the free fermion
velocity $v_{0}$ and the long-time velocity $v_{2}$ in the interacting case
converge, which is demonstrated in Fig. 4(c). For comparison, we also choose
DMS states as initial states in the non-interacting case. Again, the particles
are aligned in a periodic way (see Fig. 1(b)), with $D$ the distance between
particles within a period, in analog with $R$ in the interacting case. It is
shown that the converged value of $v_{0}$ is independent with $D$ and filling
$\rho$, as expected for free fermions. Hereafter, the velocities $v_{0}$ and
$v_{2}$ are defined according to their thermodynamic limits. We observe that
the velocity $v_{2}$ matches the dispersion velocity $v_{Q}$ extracted from
the long-time oscillations. The linear fit in Fig. 4(d) suggests the relation
$v_{0}/v_{2}=1-R\rho$. This is verified in our numerical calculation for a
broad range of $R$ and $\rho$.

The renormalization of the velocity in the interacting dynamics as compared to
the non-interacting case can be understood from the mapping between the
logical and physical basis. In this mapping, a length scale ($l$) in the
physical basis (interacting) shrinks to $l-lR\rho$ in the logical basis (non-
interacting) on average. Such shrinking has been used in the previous study on
the ground state of the 1D strongly interacting fermions for the Luttinger
parameter [24]. In our study of the dynamics, since the long-time behavior of
the strongly interacting model is captured by the plane-wave picture, it is
reasonable to expect the velocity renormalizes according to the length scale,
which then implies

$\frac{v_{0}}{v_{2}}\simeq 1-R\rho.$ (6)

This is exactly the same as what we observed in our numerical simulations.

Conclusion.— In this work, we study the quantum transport properties of
strongly interacting spinless fermions in one dimension. This model is
simulated efficiently using a fast fermion sampling algorithm at the strong
interaction limit. Despite the strong interaction, we observe robust ballistic
transport for different fillings and interaction ranges. We propose a plane-
wave description for the ballistic transport of the strongly interacting
fermions, which captures the long-time dynamics. The predominant effect of the
strong interaction on the transport is to introduce velocity renormalization.
Our results indicate the subdiffusive transport observed for the weakly
interacting spin chain in the experiment [7] would have a crossover to
ballistic behavior when approaching the strongly interacting regime. This is
worth further experimental investigation.

Acknowledgement.— We acknowledge support by National Program on Key Basic
Research Project of China (Grant No. 2021YFA1400900), National Natural Science
Foundation of China (Grants No. 11934002), Shanghai Municipal Science and
Technology Major Project (Grant No. 2019SHZDZX01), Shanghai Science Foundation
(Grants No.21QA1400500, 19ZR1471500).

## References

  * Bertini  _et al._ [2021] B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen, R. Steinigeweg, and M. Žnidarič, Finite-temperature transport in one-dimensional quantum lattice models, Rev. Mod. Phys. 93, 025003 (2021).
  * Takahashi [2005] M. Takahashi,  _Thermodynamics of one-dimensional solvable models_ (Cambridge university press, 2005).
  * Bulchandani  _et al._ [2021] V. B. Bulchandani, S. Gopalakrishnan, and E. Ilievski, Superdiffusion in spin chains, Journal of Statistical Mechanics: Theory and Experiment 2021, 084001 (2021).
  * Giamarchi [2003] T. Giamarchi,  _Quantum physics in one dimension_ , Vol. 121 (Clarendon press, 2003).
  * Hild  _et al._ [2014] S. Hild, T. Fukuhara, P. Schauß, J. Zeiher, M. Knap, E. Demler, I. Bloch, and C. Gross, Far-from-equilibrium spin transport in heisenberg quantum magnets, Phys. Rev. Lett. 113, 147205 (2014).
  * Ronzheimer  _et al._ [2013] J. P. Ronzheimer, M. Schreiber, S. Braun, S. S. Hodgman, S. Langer, I. P. McCulloch, F. Heidrich-Meisner, I. Bloch, and U. Schneider, Expansion dynamics of interacting bosons in homogeneous lattices in one and two dimensions, Phys. Rev. Lett. 110, 205301 (2013).
  * Jepsen  _et al._ [2020] P. N. Jepsen, J. Amato-Grill, I. Dimitrova, W. W. Ho, E. Demler, and W. Ketterle, Spin transport in a tunable heisenberg model realized with ultracold atoms, Nature 588, 403 (2020).
  * Eisert  _et al._ [2015] J. Eisert, M. Friesdorf, and C. Gogolin, Quantum many-body systems out of equilibrium, Nature Physics 11, 124 (2015).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239 (2016).
  * Rigol  _et al._ [2007] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Relaxation in a completely integrable many-body quantum system: An ab initio study of the dynamics of the highly excited states of 1d lattice hard-core bosons, Phys. Rev. Lett. 98, 050405 (2007).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008).
  * Narasimhan [1999] T. N. Narasimhan, Fourier’s heat conduction equation: History, influence, and connections, Reviews of Geophysics 37, 151 (1999).
  * Zotos  _et al._ [1997] X. Zotos, F. Naef, and P. Prelovsek, Transport and conservation laws, Phys. Rev. B 55, 11029 (1997).
  * Prosen [2011] T. Prosen, Exact nonequilibrium steady state of a strongly driven open xxz chain, Phys. Rev. Lett. 107, 137201 (2011).
  * Castella  _et al._ [1995] H. Castella, X. Zotos, and P. Prelovšek, Integrability and ideal conductance at finite temperatures, Phys. Rev. Lett. 74, 972 (1995).
  * Ilievski and De Nardis [2017] E. Ilievski and J. De Nardis, Microscopic origin of ideal conductivity in integrable quantum models, Phys. Rev. Lett. 119, 020602 (2017).
  * Karrasch  _et al._ [2014] C. Karrasch, J. E. Moore, and F. Heidrich-Meisner, Real-time and real-space spin and energy dynamics in one-dimensional spin-$\frac{1}{2}$ systems induced by local quantum quenches at finite temperatures, Phys. Rev. B 89, 075139 (2014).
  * De Nardis  _et al._ [2020] J. De Nardis, M. Medenjak, C. Karrasch, and E. Ilievski, Universality classes of spin transport in one-dimensional isotropic magnets: The onset of logarithmic anomalies, Phys. Rev. Lett. 124, 210605 (2020).
  * Mazur [1969] P. Mazur, Non-ergodicity of phase functions in certain systems, Physica 43, 533 (1969).
  * Scalapino  _et al._ [1993] D. J. Scalapino, S. R. White, and S. Zhang, Insulator, metal, or superconductor: The criteria, Phys. Rev. B 47, 7995 (1993).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * Fukuhara  _et al._ [2013] T. Fukuhara, P. Schauß, M. Endres, S. Hild, M. Cheneau, I. Bloch, and C. Gross, Microscopic observation of magnon bound states and their dynamics, Nature 502, 76 (2013).
  * Gómez-Santos [1993] G. Gómez-Santos, Generalized hard-core fermions in one dimension: An exactly solvable luttinger liquid, Phys. Rev. Lett. 70, 3780 (1993).
  * Li  _et al._ [2019] X. Li, G. Zhu, M. Han, and X. Wang, Quantum information scrambling through a high-complexity operator mapping, Phys. Rev. A 100, 032309 (2019).
  * Sun  _et al._ [2021] H. Sun, J. Zou, and X. Li, Fermion sampling made more efficient (2021), arXiv:2109.07358 [quant-ph] .