# Source Reconstruction for Spatio-Temporal Physical Statistical Models

Connie Okasaki
Quantitative Ecology and Resource Management Program, U of Washington,
Mevin B. Hooten
Department of Statistics and Data Sciences, The University of Texas at Austin,
and
Andrew M. Berdahl
School of Aquatic and Fisheries Sciences, U of Washington  This material is
based upon work supported by the National Science Foundation Graduate Research
Fellowship under Grant No. DGE-1762114AMB was supported by the H. Mason Keeler
Endowed Professorship in Sports Fisheries Management.

###### Abstract

In many applications, a signal is deformed by well-understood dynamics before
it can be measured. For example, when a pollutant enters a river, it
immediately begins dispersing, flowing, settling, and reacting. If the
pollutant enters at a single point, its concentration can be measured before
it enters the complex dynamics of the river system. However, in the case of a
non-point source pollutant, it is not clear how to efficiently measure its
source. One possibility is to record concentration measurements in the river,
but this signal is masked by the fluid dynamics of the river. Specifically,
concentration is governed by the advection-diffusion-reaction PDE, with an
unknown source term. We propose a method to statistically reconstruct a source
term from these PDE-deformed measurements. Our method is general and applies
to any linear PDE. This method has important applications in the study of
environmental DNA and non-point source pollution.

Keywords: Inverse Problems, GMRF, Gaussian process, Basis function, SPDE

##  1 Introduction

In many applications, a substance is emitted into a dynamic environment such
as a stream or the atmosphere, where it is then passively transported.
Measuring concentrations of such a substance is a natural monitoring tool, but
inferring emissions from concentrations requires unraveling the dynamics of
the transport medium. Often, concentrations can be modeled as the solution to
a partial differential equation (PDE), in which case the emissions can be
modeled as a source term.

We present a method using stochastic PDEs (SPDEs) to calculate the conditional
distribution of the source term of a linear PDE, given measurements of the PDE
solution. These calculations are efficient when the source can be represented
as a Gaussian Markov Random Field (GMRF). Our technique builds on the seminal
SPDE work of Lindgren et al. (2011), which suggested advection-diffusion
modeling as a potential application of their approach. Few studies have
extended this application within the SPDE context. The closest in spirit are
Sigrist et al. (2012), which used an equivalent integro-difference equation;
and Sigrist et al. (2015), which used Fourier methods to build a spatio-
temporal advection-diffusion model. However, these studies largely demonstrate
the case of a single drift vector rather than an entire advection field.
Outside the SPDE context, spatial and spatio-temporal statistical models of
concentration within dynamic mediums are more common, but often either ignore
transport and emissions effects (e.g. Moraga et al. 2017), or incorporate
these effects in a non-mechanistic (e.g. Cameletti et al. 2019) or pseudo-
mechanistic manner (e.g. Ver Hoef et al. 2006). Here we use explicit
advection-diffusion mechanisms, within the SPDE context in order to obtain
sparse precision matrices, and can incorporate a spatio-temporally varying
advection field.

Our method provides a mechanistic kernel for modeling concentrations, as well
as a way to infer emissions from those concentrations. The latter is novel in
the field of spatial statistics, and statistically more difficult, because
transport distorts the original emissions signal. We therefore focus on the
latter application. In Section 2, we provide background information: a review
of the SPDE approach for generating the Matérn covariance, a review of the
finite element method as it is used in this approach, and an overview of the
advection-diffusion-reaction equation. In Section 3, we outline the method we
used for source reconstruction, as well as the method we used for quantifying
the accuracy of our source reconstructions, a method for constructing spatio-
temporal versions of our model, the MCMC method we used for sampling, and
finally any specifics regarding our three case studies. In Section 4, we
describe the results of these three case studies.

###  1.1 Motivation

Our focus on inferring emissions is motivated by several continuing and
emerging challenges in environmental science. The first and most direct
application is the problem of non-point source air and water pollution. There
is a long history of investigation into the inverse source problem in these
domains (e.g. Panofsky 1969), but study has been largely focused on the more
mathematically tractable point-source case (e.g., Stockie 2011). Nevertheless,
non-point sources are a major concern for regulators, in part because the
traditional policies used to regulate pollution are difficult to implement for
these sources (Xepapadeas, 2011).

Another major motivation of this work is environmental DNA111Environmental DNA
is particles of DNA shed from an organism into the environment. These
particles can be used to passively detect the presence of particular species
(Rees et al., 2014), or in certain circumstances can be used to estimate
biomass (Rourke et al., 2022). In theory, provided accurate parameters
describing the physical-chemical dynamics of environmental DNA, one could use
it to derive a spatio-temporal map of biomass using only non-invasive
sampling. (eDNA). The eDNA literature has made great strides in the last
decade and techniques such as DNA metabarcoding (Deiner et al., 2017) have
seen wide acceptance. More recently, scientists have investigated the use of
quantitative polymerase chain reaction analysis (qPCR222The polymerase chain
reaction is a biological reaction used to amplify a particular DNA sequence if
it is present in a sample. It is traditionally used as a detector of
presence/absence. However, in recent years qPCR has been developed as a
technique to estimate not just presence/absence but the _amount_ of a
particular DNA sequence present in the original sample.) to estimate the
abundance of organisms (e.g., Doi et al. 2015,Doi et al. 2017). Once the
dynamic properties of eDNA (shedding rate, decay rate, etc) can be quantified
in field settings, eDNA offers a promising tool for non-invasive abundance
estimation (Harrison et al., 2019). Our method offers a mathematical framework
in which many different dynamic properties can be incorporated to infer
spatial and spatio-temporal abundance.

The inverse source problem is broad, and includes many problems beyond these
two motivating applications. Other literature in this field includes inferring
sources of heat (El Badia & Ha-Duong, 2002); acoustic waves (Alves et al.,
2009); and electromagnetic signals, particularly EEG data (Grech et al.,
2008). We do not consider these applications in this paper, but the adaptation
of our method to these domains is a promising area of future study.

##  2 Background

###  2.1 Statistical Background

We assume that the source, or emissions, $f(\bm{{s}})\sim
GP(\mu(\bm{{s}}),k(\bm{{s}},\bm{{s}}^{\prime}))$ is a Gaussian process with
mean function $\mu$ and covariance function $k$. We assume that $f(\bm{{s}})$
cannot be measured directly, and is the primary target of inference.
Furthermore, we assume that $f$ can be represented using a Gaussian Markov
Random Field (GMRF; for general reference see Rue & Held (2005)) so that for
any particular finite set of points $\bm{{s}}$,
$\bm{{f}}(\bm{{s}})\sim\mathcal{N}(\bm{{\mu}}(\bm{{s}}),\bm{{Q}}(\bm{{s}},\bm{{s}})^{-1})$
with a sparse precision matrix $\bm{{Q}}$. In the special case that $k$ is the
Matérn covariance function,

$k(\bm{{s}},\bm{{s}}^{\prime})=\frac{\sigma^{2}}{2^{\nu-1}\Gamma(\nu)}(\kappa||\bm{{s}}^{\prime}-\bm{{s}}||)^{\nu}K_{\nu}(\kappa||\bm{{s}}^{\prime}-\bm{{s}}||),$
(1)

Lindgren et al. (2011) provided a method for directly generating a sparse
$\bm{{Q}}$ for certain integer and half-integer values of the smoothness
parameter $\nu$. Their approach for doing so was based on applying the finite
element method (FEM) to a particular set of nested diffusion equations whose
solutions are governed by the Matérn covariance function. We use the same
method to translate from emissions to concentration, using a different PDE
(see Equation 2).

###  2.2 The Advection-Diffusion-Reaction Equation

We now consider the advection-diffusion-reaction equation in more detail.
Suppose that we seek to infer the source of a chemical, but only have
measurements of chemical concentration in a fluid medium. The chemical flows
with, and disperses through, that medium, and possibly decays linearly over
time. The chemical concentration is then described by the equation:

$\dfrac{\partial u(\bm{{s}},t)}{\partial
t}+\nabla\cdot(\bm{{v}}(\bm{{s}},t)u(\bm{{s}},t))-\nabla\cdot(D\nabla
u(\bm{{s}},t))+ru(\bm{{s}},t)=f(\bm{{s}},t),$ (2)

where $\nabla$ is the del operator for the spatial dimensions, $D$ is the
diffusion coefficient, $r$ is the decay rate, $\bm{{v}}(\bm{{s}},t)$ is the
advection field, assumed here to be known, $f(\bm{{s}},t)$ is the source
function, a general Gaussian process in space and time, and $u(\bm{{s}},t)$ is
the solution, or concentration function. If we assume that
$f(\bm{{s}},t)=f(\bm{{s}})$ is homogeneous in time and that
$u(\bm{{s}},t)=u(\bm{{s}})$ has reached steady state, the PDE simplifies to

$\nabla\cdot(\bm{{v}}(\bm{{s}})u(\bm{{s}}))-\nabla\cdot(D\nabla
u(\bm{{s}}))+ru(\bm{{s}})=f(\bm{{s}}).$ (3)

For example, this model may be appropriate when trying to estimate the origin
of non-point source pollution (e.g., runoff pollution); the source of isotopes
found in water samples333Many elements found in nature have multiple naturally
occurring and scientifically useful isotopes, such as tritium (3H),
radiocarbon (14C), and oxygen-18 (18O). Such isotopes can be naturally
enriched by environmental processes such as evaporation and can therefore be
used to identify the source of water flow (groundwater vs. snow pack vs.
rainfall; e.g. Hao et al. 2019) or match a sample with a particular location
(e.g. Rachel et al. 2008).; or the biomass of fish from eDNA measurements.
This inverse problem has been considered in applied math studies, and is
usually solved either in the case where the chemical is assumed to originate
with unknown strength from a small finite number of unknown point locations
(El Badia et al., 2005; El Badia & Hamdi, 2007; Lushi & Stockie, 2010), or
using regularization techniques from the inverse problems literature that do
not admit a probabilistic interpretation and do not include an estimate of
uncertainty (Porter & Devaney, 1982; Engl et al., 1996; Yan et al., 2008).

Previous statistical work has examined the advection-diffusion equation in
similar contexts using integro-difference equations (Sigrist et al., 2012) and
Fourier methods (Sigrist et al., 2015). Stroud et al. (2010) used the finite
difference method to account for these practicalities but used an ensemble
Kalman filter method to conduct inference, rather than approaching the problem
from an SPDE perspective.

###  2.3 The Finite Element Method

In general, discretization methods translate partial differential equations
into systems of linear equations

$\bm{{K}}\bm{{u}}=\bm{{L}}\bm{{f}}.$ (4)

In the finite element method, $\bm{{K}}$ is known as the stiffness matrix and
$\bm{{L}}$ is known as the mass matrix. In the finite difference method, often
$\bm{{L}}=\bm{{I}}$. Essential boundary conditions (i.e., those that are not
automatically satisfied; usually Dirichlet) modify this basic equation. We
demonstrate this modification in the Appendix. We refer to $\bm{{u}}$ as the
solution vector and $\bm{{f}}$ as the source vector, with precision matrices
$\bm{{Q}}_{u}$ and $\bm{{Q}}_{f}$.

Our methodology for generating a sparse precision matrix $\bm{{Q}}_{u}$ given
$\bm{{Q}}_{f}$ is essentially the same as in Lindgren et al. (2011). We
generalize their notation and define $\mathcal{L}$ to be any linear partial
differential operator and $\mathcal{B}$ to be a general Dirichlet, Neumann, or
Robin boundary condition, to produce a general formulation of a PDE

$\begin{split}\mathcal{L}u(\bm{s})&=f(\bm{s})\qquad\mbox{for
}\bm{{s}}\in\Omega\\\ \mathcal{B}u(\bm{s})&=c(\bm{s})\qquad\mbox{for
}\bm{{s}}\in\partial\Omega.\end{split}$ (5)

Neumann boundary conditions are the natural boundary conditions in our
analyses and we use them throughout, but Dirichlet boundary conditions could
also be used, as could non-zero Neumann boundary conditions and Robin boundary
conditions. Boundary conditions may also be specified in terms of lower-
dimensional random fields (see Appendix for details).

The finite element method discretizes Equation 5 by finding a weak solution
using a specified set of test and trial functions. In Lindgren et al. (2011),
and in this study, both classes are spanned by a series of piecewise linear
“bump” functions defined to be 0 over most of the domain, and 1 at one
specific node chosen from a set of finite points in space. Other classes may
be advantageous for particular SPDEs, including the advection-diffusion
equation (which in future work may benefit from the Streamline-Upwind Petrov
Galerkin bases; Brooks & Hughes 1982) and Maxwell’s equations (which benefit
from the use of specialized bases such as the Nédélec elements; Monk et al.
2003). For a more comprehensive review of the FEM in general see Brenner et
al. (2008), or for an accessible review of the FEM applied to the Matérn
equation see Bakka (2018).

Conceptually, the finite element method serves to translate the PDE in
Equation 5 into a matrix equation such as Equation 4. Classically, one then
solves for $\bm{{u}}=\bm{{K}}^{-1}\bm{{L}}\bm{{f}}$. Lindgren et al. (2011)
observed that in doing so, we imply that if $\bm{{f}}$ is multivariate normal
(MVN) with mean $\bm{{\mu}}_{f}$ and precision $\bm{{Q}}_{f}$ then $\bm{{u}}$
is MVN with mean $\bm{{K}}^{-1}\bm{{L}}\bm{{\mu}}_{f}$ and precision
$\bm{{K}}^{T}\bm{{L}}^{-T}\bm{{Q}}_{f}\bm{{L}}^{-1}\bm{{K}}$. To ensure
sparsity of the resulting matrix, we replace $\bm{{L}}$ with a diagonal matrix
$\tilde{\bm{{L}}},$ formed by summing the rows of $\bm{{L}}$ and placing the
resulting sums along the diagonal.

##  3 Methods

###  3.1 Finite Element Calculations

Although the finite difference method would produce conceptually identical
results, we use the finite element method to compute discretized spatial
derivatives. Because the stiffness and mass matrices are standard components
of finite element analysis in engineering, they can be calculated by
leveraging existing FEM software such as the Finite Element Computational
Software (FEniCS (Alnæs et al., 2015). We rely on FEM software to implement
our stiffness and mass matrices, to improve reproducibility and facilitate
generalizability. We assume that the FEM approach is appropriate for any SPDE
in which the stochastic source term is well-behaved and in which the FEM
equation converges to the true solution in the deterministic case.

###  3.2 Source Reconstruction

Our basic observation is that physical meaning can sometimes be ascribed to
$f(\bm{{s}})$ and that the finite element representation $\bm{{f}}$ can be
calculated by inverting Equation 4 in the reverse direction
$\bm{{f}}=\bm{{L}}^{-1}\bm{{K}}\bm{{u}}$. Our method for source reconstruction
given a set of parameters (which generally must be estimated) therefore
follows essentially the same logic as classical kriging. We assume that the
solution $\bm{{u}}$ is observed through some set of linear functionals which
can be discretized to a matrix $\bm{{A}}$. We therefore observe
$\bm{{y}}=\bm{{Au}}+\bm{{\epsilon}}$ with i.i.d. Gaussian noise. The
conditional field $\bm{{u}}|\bm{{y}}$ is therefore

$\bm{{u}}|\bm{{y}}\sim
N\left(\bm{{\mu}}_{u}+\frac{1}{\sigma^{2}}(\bm{{Q}}_{u}+\frac{1}{\sigma^{2}_{\epsilon}}\bm{{A}}^{T}\bm{{A}})^{-1}\bm{{A}}^{T}(\bm{{y}}-\bm{{A}}\bm{{\mu}}_{u}),(\bm{{Q}}_{u}+\frac{1}{\sigma_{\epsilon}^{2}}\bm{{A}}^{T}\bm{{A}})^{-1}\right).$
(6)

This is a standard kriging computation for a GMRF (Eq. 2.15 of Rue & Held
2005). However, we propose to take the additional step of kriging
$\bm{{f}}|\bm{{y}}$ by observing that
$\bm{{f}}|\bm{{y}}=\bm{{L}}^{-1}\bm{{K}}\bm{{u}}|\bm{{y}}$ and therefore

$\bm{{f}}|\bm{{y}}\sim
N\left(\bm{{L}}^{-1}\bm{{K}}\bm{{\mu}}_{u|y},\bm{{L}}^{-1}\bm{{K}}(\bm{{Q}}_{u}+\frac{1}{\sigma_{\epsilon}^{2}}\bm{{A}}^{T}\bm{{A}})^{-1}\bm{{K}}^{T}\bm{{L}}^{-1}\right).$
(7)

To facilitate numerical calculation we use the Cholesky factors of
$\bm{{Q}}_{u}$ and
$\bm{{Q}}_{u}+\frac{1}{\sigma_{\epsilon}^{2}}\bm{{A}}^{T}\bm{{A}}$. We
calculate these Cholesky factors using CHOLMOD functions (Chen et al., 2008)
provided by the sksparse library in Python. Specifically, we calculate the
Cholesky factor for $\bm{R}^{T}\bm{R}=\bm{{Q}}_{u}$, and use the sparse
Cholesky update methods provided by CHOLMOD to update that factor to account
for the addition of a low-rank matrix
$\frac{1}{\sigma^{2}_{\epsilon}}\bm{{A}}^{T}\bm{{A}}$. Note that this update
step is computationally difficult when implementing our model in space-time.

Linear regression on $\bm{{f}}$ can be accommodated within this kriging
framework; this calculation is shown in the Appendix.

###  3.3 Accuracy Quantification

To determine the accuracy of our method for mechanistic kriging, we quantified
error using the $L^{2}$ norm over the domain (excluding any buffer regions)
relative to the kriging estimator:

$L^{2}(x|y)=\sqrt{\int_{\Omega}(x|y-\mu_{x|y})^{2}d\Omega}.$ (8)

In simulations, we calculated the integral over the finite element basis
functions using FEniCS native integration functions. However, assuming the
mesh has cells of approximately equal extent, and assuming that all
hyperparameters are known perfectly so that all error is produced by Gaussian
process variability and observation error, we may approximate the error as

$L^{2}(x|y)\approx\sqrt{\frac{V}{M}\rm{tr}\left(\bm{{I}}_{\rm
int}\bm{{\Sigma}}_{\bm{{x}}|\bm{{y}}}\bm{{I}}_{\rm int}\right)},$ (9)

where $\bm{{x}}$ represents either $\bm{{f}}$ or $\bm{{u}}$, $V$ is the volume
of the interior region $\Omega$, $\bm{{\Sigma}}_{\bm{{x}}|\bm{{y}}}$ is the
conditional covariance of vector $\bm{{x}}$ given the vector $\bm{{y}}$ of
observations of $\bm{{u}}$, $\bm{{I}}_{\rm int}$ is a diagonal matrix whose
entries are the fractions of each basis function contained within $\Omega$,
and $M$ is $\rm{tr}(\bm{{I}}_{\rm int}).$

Additionally, to quantify the rate of convergence of error, we consider the
possibility of observing the same set of of locations $\bm{{A}}^{T}\bm{{A}}$ a
total of $n$ times, resulting in a precision matrix

$\bm{{Q}}_{\bm{{u}}}+\frac{n}{\sigma^{2}_{\epsilon}}\bm{{A}}^{T}\bm{{A}}.$

We relax the interpretation that we actually observe locations an integral
number of times so as to allow a derivative with respect to $n$. We
reparameterize (with $m$ the number of rows in $\bm{{A}}$), and set
$\zeta=\log(mn)$; we also reparameterize to obtain the log of the error.
Differentiating the log-error with respect to $\zeta$, we can obtain the local
polynomial convergence rate. Further, if we do not use a buffer zone, thereby
allowing us to observe all nodes, we can also show that as $n\to\infty$, the
polynomial convergence rate approaches $-1/2$ for both the source and solution
functions. This corresponds to a $1/\sqrt{N}$ type asymptotic convergence for
the $L^{2}$ error. We show that this asymptotic behavior is not observed
empirically, either due to the inclusion of a buffer region, or due to slow
convergence to the asymptote. Instead, the local polynomial convergence rate
calculation is more useful in practice. Full calculations are shown in the
Appendix.

###  3.4 Spatio-Temporal Modeling

Our implementation of time-dependence relies on a finite element method in
space and a finite difference method in time. This results in a sparse vector
auto-regressive (VAR) model for the spatial variables (Wikle & Hooten, 2010).
We discretize the time interval $[t_{0},t_{1}]$ into small units of $\Delta
t=(t_{1}-t_{0})/N$. We can then derive a matrix equation from any arbitrary
SPDE

$\dfrac{\partial u(\bm{{s}},t)}{\partial
t}+\mathcal{L}(t)u(\bm{{s}},t)=f(\bm{{s}},t)$ (10)

(where $\mathcal{L}(t)$ is an arbitrary partial differential operator; for
example, in Equation 2, we have
$\mathcal{L}(t)u=\nabla\cdot(\bm{{v}}u)-\nabla\cdot(D\nabla u)+ru$) by first
applying a backward Euler discretization in time

$\begin{split}\dfrac{u(\bm{{s}},t+\Delta t)-u(\bm{{s}},t)}{\Delta
t}+\mathcal{L}u(\bm{{s}},t+\Delta t)&\approx f(\bm{{s}},t+\Delta t),\\\
u(\bm{{s}},t+\Delta t)+\Delta t\mathcal{L}u(\bm{{s}},t+\Delta t)&\approx
u(\bm{{s}},t)+\Delta tf(\bm{{s}},t+\Delta t),\end{split}$ (11)

and then applying an FEM discretization in space

$\begin{split}\bm{{L}}\bm{{u}}_{t+\Delta t}+\Delta t\bm{{K}}_{t+\Delta
t}\bm{{u}}_{t+\Delta t}&\approx\bm{{L}}\bm{{u}}_{t}+\Delta
t\bm{{L}}\bm{{f}}_{t+\Delta t},\\\ (\bm{{I}}+\Delta
t\bm{{L}}^{-1}\bm{{K}}_{t+\Delta t})\bm{{u}}_{t+\Delta
t}&\approx\bm{{u}}_{t}+\Delta t\bm{{f}}_{t+\Delta t},\\\
u(\bm{{s}},t)&\approx(\bm{{I}}+\Delta t\bm{{L}}^{-1}\bm{{K}}_{t+\Delta
t})^{-1}(\bm{{u}}_{t}+\Delta t\bm{{f}}_{t+\Delta t}),\end{split}$ (12)

where all functions of space and time are replaced by vectors in space,
subscripted by time, and the operator $\mathcal{L}(t+\Delta t)$ is discretized
into $\bm{{K}}_{t+\Delta t}$. This equation can be solved inductively by

$\begin{bmatrix}\bm{{u}}_{1}\\\ \bm{{u}}_{2}\\\ \vdots\\\
\bm{{u}}_{N}\end{bmatrix}=\begin{bmatrix}\bm{{M}}_{1}\\\
\bm{{M}}_{2}\bm{{M}}_{1}\\\ \vdots\\\
\bm{{M}}_{N}...\bm{{M}}_{1}\end{bmatrix}\bm{{u}}_{0}+\Delta
t\left[\begin{array}[]{cccc}\bm{{M}}_{1}&&&\\\
\bm{{M}}_{2}\bm{{M}}_{1}&\bm{{M}}_{2}&&\\\ \vdots&\vdots&\ddots&\\\
\bm{{M}}_{N}...\bm{{M}}_{1}&\bm{{M}}_{N}...\bm{{M}}_{2}&\cdots&\bm{{M}}_{N}\\\
\end{array}\right]\begin{bmatrix}\bm{{f}}_{1}\\\ \bm{{f}}_{2}\\\ \vdots\\\
\bm{{f}}_{N}\end{bmatrix},$ (13)

where $\bm{{M}}_{i}=(\bm{{I}}+\Delta t\bm{{L}}^{-1}\bm{{K}}_{t_{0}+i\Delta
t})^{-1}$. The above Toeplitz matrix has sparse inverse

$\bm{{R}}=\frac{1}{\Delta t}\left[\begin{array}[]{cccc}\bm{{M}}_{1}&&&\\\
-\bm{{I}}&\bm{{M}}_{2}&&\\\ &\ddots&\ddots&\\\ &&-\bm{{I}}&\bm{{M}}_{N}\\\
\end{array}\right].$ (14)

A forward Euler discretization yields a very similar matrix, but requires that
the time step $\Delta t$ be small enough to ensure numerical stability.

Suppose now that the full, discretized, spatio-temporal process $\bm{{f}}$ has
sparse precision matrix $\bm{{Q}}_{f}$ and that we have a deterministic
initial condition $\bm{{u}}_{0}=\bm{{0}}$ (accounted for by including a
temporal buffer region before the region of interest). Then the precision
matrix for $\bm{{u}}$ is $\bm{{Q}}_{u}=\bm{{R}}^{T}\bm{{Q}}_{f}\bm{{R}}$. Our
initial condition requires us to extend the temporal discretization prior to
the first time-of-interest, to avoid boundary effects. Alternatively, the
steady-state equation could be used to derive an approximate spatial prior for
the initial condition, although this should still be coupled with a (perhaps
smaller) buffer region.

This model uses sparse precision matrices for efficient computation. The block
bi-diagonal Cholesky factors allow rapid computation in a high dimensional
space, in part because the blocks are themselves sparse matrices. This makes
computing likelihoods efficient, scaling linearly in the number of time-steps
modeled. The computational difficulty of this model arises from the Cholesky
update step (i.e., the factorization of
$\bm{Q}_{\bm{u}}+\frac{1}{\sigma^{2}_{\epsilon}}\bm{A}^{T}\bm{A}$), necessary
both for inferring the distribution for the emissions or concentration maps,
conditional on the observations (the kriging step) and for calculating the
marginal likelihood function (the marginalization step). These conditional
updates break the structure of our Cholesky factors, resulting in much greater
computational burden. In essence, the problem is that information from each
observation is transported over time, resulting in long-range dependencies. In
densely sampled systems, this problem might be solved by calculating the
temporal range of the covariance function, and truncating dependencies after a
fixed number of time-steps, thus limiting the bandwidth of the resulting
updated Cholesky factor.

####  3.4.1 Nested Diffusion Definition of Spatio-temporal Matérn Source

The need for a non-mechanistic spatio-temporal source distribution compatible
with our method motivates the following definition of a three-parameter
spatio-temporal Matérn field

$\left(\tau\dfrac{\partial}{\partial
t}+\kappa^{2}-\Delta\right)^{\alpha/2}u(\bm{{s}},t)=\mathcal{W}(\bm{{s}},t).$
(15)

This is a nested spatio-temporal diffusion equation, where the more general
source term $f(\bm{{s}},t)$ is specified to be constant-in-time spatial white
noise (denoted $\mathcal{W}(x)$ for consistency with Lindgren et al. 2011)
then the spatial Matérn SPDE is the steady state solution. Although the Matérn
model does not necessarily require interpretation, the interpretation of this
nested diffusion model is more straightforward than an alternative model for
which the spatial Matérn SPDE is also the steady state solution

$\tau\dfrac{\partial u(\bm{{s}},t)}{\partial
t}+\left(\kappa^{2}-\Delta\right)^{\alpha/2}u(\bm{{s}},t)=\mathcal{W}(\bm{{s}},t).$
(16)

The two equations are equivalent when $\alpha=2$; both can be seen as a
diffusion-decay equation with natural time scale $\tau$. In this circumstance,
“randomness” in the form of white noise is spread out across space by means of
a familiar mechanism. At higher values of $\alpha$ the equations and their
interpretations diverge. Equation 15 can be seen to be a nested form of the
diffusion-decay equation; converting from an Eulerian to a Lagrangian
perspective, one can imagine random particles entering a system according to
$\mathcal{W}(\bm{{s}},t)$ and diffusing, only to continue emitting secondary
particles that themselves diffuse, and so on. In contrast Equation 16, even
for integer values of $\alpha/2$, involves at minimum fourth-order-in-space
operators coupled to a first-order-in-time derivative. Fourth-order equations,
while not at all uncommon in specialized mathematical fields, are relatively
exotic. This definition aligns with that proposed in (Bakka et al., 2020),
with $\alpha_{s}=2$ and $\mathcal{E}_{Q}$ white in both space and time.

If we wish to avoid fractional diffusion we must assume $\alpha=2n$ for
$n\in\mathbb{Z}_{+}$, for which we obtain

$\bm{{Q}}_{u}^{(2n)}=\frac{\bm{{R}}^{T}}{\Delta
t}^{n}\left(\bm{{L}}\otimes\bm{{I}}\right)\frac{\bm{{R}}}{\Delta t}.$ (17)

This results in $\bm{{Q}}_{f}=\bm{{L}}\otimes\bm{{I}},$ which is
straightforward to calculate. For example, Figure 1 shows the results of a
numerical simulation of this distribution. Although a spatial buffer-zone was
included in these calculations and removed from the plot, the temporal buffer-
zone was not removed, to demonstrate the transient behavior imposed by our
deterministic initial condition. A small zone with length on the order $\tau$
displays lower spatio-temporal variance than the remainder of the plot, but no
differences appear to persist beyond that. For this set of parameters, we used
$\Delta t=0.05$ and $2000$ time steps. With 751 spatial nodes, simulation
required $\approx 21$ seconds on a desktop computer with a 3.2 GHz 6-core
processor and 32 GB of RAM.

Figure 1: The results of a numerical simulation of the proposed spatio-
temporal Matérn distribution. Note that this function is both the solution to
an SPDE in one context ($u(s,t)$ in Equation 15) but that we propose to use it
as a source function in another context ($f(s,t)$ in Equation 10). Parameters
used were: $\tau=2,\alpha=4,\kappa=1$, corresponding to $\rho\approx 7$. A
spatial buffer zone of length 15 on each side was excluded from this figure.
No temporal buffer zone was removed. Space-time was discretized using 751
spatial nodes, and $10^{4}$ time steps.

###  3.5 Parameter Estimation

To estimate the parameters of our air pollution model (see Section 3.6.3) we
used MCMC. To facilitate a Gibbs update for the variance, we specified
inverse-gamma priors for $\sigma^{2}_{f}$ the variance of the emissions
function, and we specified gamma priors for $\rho,v,D,$ and $r$ (respectively,
the range; ratio of observation variance $\sigma^{2}_{\epsilon}$ to emissions
variance $\sigma^{2}_{f}$; diffusion coefficient; and decay coefficient). All
other parameters are updated in Metropolis-Hastings steps.

We marginalized over the concentration and emission vectors $\bm{{u}}$ and
$\bm{{f}}$, but used a second-stage sampler to draw from the full-conditional
distributions of $\bm{{u}}$ and $\bm{{f}}$ to compute their posterior means
and pointwise standard deviations.

To treat the range as a parameter we used the empirical equation
$\rho=\sqrt{(8\nu)}/\kappa$ which relates range to the Matérn parameters
(Lindgren et al., 2011) and transformed the log-likelihoods accordingly. To
use Gibbs sampling for $\sigma^{2}_{f}$ we reparameterized to the variance
ratio: $v=\frac{\sigma^{2}_{\epsilon}}{\sigma^{2}_{f}}$. We transformed the
log-likelihoods accordingly. Due to the identifiability problems with the
Matérn $\alpha$ parameter (Lindgren et al., 2011), we assumed a priori that
$\alpha=2$ and therefore $\nu=1$, and that the prior mean of $\bm{f}$ was
everywhere zero.

To verify the performance of this MCMC algorithm, we also conducted MCMC in
our 1-D non-temporal model (see Section 3.6.1). In this model we also
reparameterized the hyperparameter
$\sigma^{2}_{\beta}=\sigma^{2}_{f}v_{\beta}$ to allow a Gibbs update for the
emissions variance. Our air pollution model did not include covariates and
therefore did not require this parameter. We also used different numbers of
iterations in the two examples. Otherwise the algorithms were identical. Note
that we did not conduct MCMC in our 1-D spatio-temporal model (see Section
4.2) due to computational constraints. In this case study, we only demonstrate
the mechanics of the kriging step.

###  3.6 Case Studies

####  3.6.1 Case Study 3.6: Inferring a 1-Dimensional Spatial Source

We conducted an analysis in one-dimension of space, for a system at
equilibrium, with a source governed by a linear model.

Suppose we make measurements of runoff pollution $u(s)$ at various points in a
1-D stream. The land alongside the stream may be used for several different
purposes; for example, it may be agricultural, residential, or industrial. The
primary source of variation in pollutant source $f(s)$ may be due to
variations in land use (represented by a factor variable), so that
$f(s)=\bm{{x}}(s)^{T}\bm{{\beta}}+\epsilon(s)$ where $\epsilon$ can be
represented as a mean-zero GMRF and $\bm{{\beta}}$ has a normal prior. In this
case study, to demonstrate the kriging step, we assumed all dynamic parameters
were known and just calculated the kriging step. We applied the finite element
method to Equation 3 to obtain $\bm{{K}}$ and $\bm{{L}}$. Then, we calculated
the joint conditional distribution for $\bm{{f}}$ and $\bm{{\beta}}$ given
measurements of $\bm{{y}}=\bm{{Au}}+\bm{{\epsilon}}$ as described in Section
3.2 by first calculating the joint conditional distribution of $\bm{{u}}$ and
$\bm{{\beta}}$ then inverting Equation 4 to obtain the conditional
distribution for $\bm{{f}}$.

To demonstrate the “kriging step” we first made predictions assuming the
dynamics were known. Then, to demonstrate the accuracy of this kriging
calculation, we calculated the average $L^{2}$ error across 30 simulations for
each of 30 sampling densities ranging from a single data point up to 200,000
data points. For each sampling density, samples were evenly distributed across
the interior of the fixed domain, but not into the buffer region. For the
purposes of approximating the error we assumed that we observed nodes
uniformly $\bm{{A}}^{T}\bm{{A}}=\bm{{I}}_{\rm int}$ a fractional number of
times equal to $n=N/m$, the true sample size divided by the number of rows of
$\bm{{A}}$. Note that the relative error that we estimate here is purely
interpolation error – all parameters are assumed known.

Finally, we conducted an MCMC analysis (described in Section 3.5) to
demonstrate the efficacy of our inference in the air pollution case study. In
this case study, we ran 4 chains for 5,000 iterations each, with a thinning
rate of 5 for a total of 4,000 samples. We also used 5,000 iterations of burn-
in. To match our air pollution case study, we specified weakly informative
priors. Since we marginalized out the source, solution, and regression
coefficients, we calculated these posterior distributions as generated
quantities after completing the main MCMC sampling.

The true parameter values used for simulation in this case study can be found
in Table 1.

Parameter | $\alpha$ | $\rho$ | $D$ | $r$ | $\sigma^{2}_{f}$ | $\sigma^{2}_{\epsilon}$ | $\sigma^{2}_{\beta}$
---|---|---|---|---|---|---|---
Case Study 3.6 | 2 | 2 | 0.75 | 0.2 | 10 | 10 | 25
Case Study 1 | 4 | 5.3 | 0.25 | 0.05 | 10 | 10 | NA
Table 1: The true parameters used for simulation in Case Studies 3.6 and 1.

####  3.6.2 Case Study 1: Inferring a 1-Dimensional Spatio-Temporal Source

In the previous section, we assumed a PDE in steady-state. Steady-state PDEs
assume that dynamics and emissions are constant, or change on a time scale
much longer than the time scale over which the steady-state solution is
reached. However, in many applications these assumptions do not hold. In these
cases, a full spatio-temporal model may be necessary. We provide a
demonstration of this model, but we note that spatio-temporal models are
computationally much more challenging, and we defer a more intensive spatio-
temporal application for future work. In particular we demonstrate only the
kriging step that is the foundation of these calculations, and not a more
intensive MCMC analysis.

In this case study, we simulate a 1-D domain, with small enough dimensionality
that computation remains feasible (101 spatial nodes over 2000 time steps,
with simulation taking $\approx$1 second on a desktop computer with a 3.4 GHz
6-core processor and 32 GB of RAM. Full spatio-temporal emissions estimation
in higher-dimensional spaces is a promising area of future study.

The true parameter values used for simulation in this case study can be found
in Table 1.

####  3.6.3 Case Study 3.6.2: Inferring Sources of Air Pollution in the U.S.

Finally, we implemented the steady-state advection-diffusion-reaction equation
to conduct source inference on PM2.5 over the United States for January 1st
2019. We obtained PM2.5 data from the EPA (US Environmental Protection Agency,
2021) and parameterized our wind-speed field using NOMADS (Rutledge et al.,
2006). We assumed a homogeneous prior mean for $\bm{{f}}$, although population
density, automobile density, land use, and many other covariates could also be
included. We included a buffer region of approximately 4 degrees on all sides,
and projected locations from degrees latitude/longitude to kilometers,
preserving the distance between longitudes at all latitudes. We also projected
advection velocities to ensure that the travel time between two points was
preserved, although distances between two points and advection speeds were in
general not preserved.

We assumed that the emissions (source function; $f(\bm{{s}})$) were governed
by a Matérn Gaussian process and that the variance and range of this process
were unknown, along with the diffusion and decay coefficients. We fit this
model using an MCMC algorithm described in Section 3.5. In this case study we
ran 4 chains for 25,000 iterations each, with a thinning rate of 5 for a total
of 20,000 samples. We also used 25,000 iterations of burn-in.

We specified weakly informative priors based on existing literature ($\rho$
(Song et al., 2018); $D$ (Byun & Schere, 2006); $r$ (Watson et al., 2000; Cao
et al., 2013); $\sigma^{2}_{f}$ (Tucker, 2000); $v$ (Peters et al., 2001)) and
understanding of the mechanisms involved. Prior parameters and quantiles, and
posterior statistics can be found in Table 2.

##  4 Results

###  4.1 1-Dimensional Spatial Source Results

Figure 2 shows the results of our kriging analysis, with three land-use zones
applying differing constant contributions to the source. This demonstrates the
ability to incorporate a linear regression model into our method and to obtain
analytical conditional distributions when dynamics are perfectly known.

Figure 2: The results of a numerical experiment for conducting source
reconstruction on the advection-diffusion-reaction equation, with a linear
regression on spatial factors. Panel (a) shows the spatially-variable flow
rate (always to the right). Panel (b) shows inference on $\bm{{u}}$. Panel (c)
shows inference on $\bm{{f}}=\bm{{X}}\bm{{\beta}}+\bm{{\epsilon}}$. Note the
discontinuities in $\bm{{f}}$ generated by discontinuities in $\bm{{X}}$. Dots
show positions of measurements and their observed values.

When quantifying the error incurred by kriging in this case study, we found
that the local polynomial convergence rate approximations for the $L^{2}$
error agreed with our empirical calculations at all but the lowest sample
sizes ($\lesssim$ 10). We also found that across a wide range of moderate
sample sizes, the error converges at a rate slower than $1/\sqrt{N}$.
Specifically, we observed a convergence rate of approximately $N^{-0.4}$ when
estimating concentrations (i.e., PDE solutions), and an error convergence rate
of approximately $N^{-0.1}$ when estimating emissions (i.e., PDE sources).
Hypothesizing that lower rates of convergence are the result of information
being obscured by the dynamics, particular by diffusion and decay, and of
observation error being too large to effectively “see past” these obfuscation,
we performed the same set of simulations with the diffusion rate ($D$), decay
rate ($\lambda$), and observation variance ($\sigma^{2}_{\epsilon}$) all
reduced by a factor of 10. In this case, we observed faster convergence rates
of $N^{-0.41}$ and $N^{-0.19}$ respectively.

|
---|---
Figure 3: Empirical vs. approximate plots showing the convergence of error for
our method when estimating both the source and solution of a PDE from
solution-only measurements in Case Study 3.6. Note that our analytical
approximations underestimate the error initially, but become accurate at a
sample size of $N<10$. In the first simulation (a), $D=0.75,r=0.2,$ and
$\sigma^{2}_{\epsilon}=5$, while in the second (b) $D=0.075,r=0.02,$ and
$\sigma^{2}_{\epsilon}=0.5$. Figure 4: The results of Case Study 3.6’s MCMC
analysis. Prior PDF had been scaled up by the number of samples to be
comparable to the histogram.  Figure 5: The results of Case Study 3.6’s MCMC
analysis. Black dotted lines represent the true simulated source and solution
functions. Red dotted lines represent the pointwise posterior mean. Background
lines represent the kriging estimate from each MCMC sample. Note that the
solution (or concentration) is estimated precisely but the source (or
emissions) is more uncertain. Both functions lie largely within the range of
variation among samples.  Figure 6: The generated distributions for the
source, solution, and regression coefficients given the posterior samples of
the parameters. The vertical line represents the true value of the
coefficient. Note that all three coefficients are underestimated on average,
perhaps due to oversmoothing in the kriging step.

Upon conducting MCMC we found that some parameters could be inferred well
while others suffered from issues of identifiability. The posteriors for the
diffusion, decay, range, and emissions variance (shown in Figure 4) improved
noticeably over the priors, although the decay showed a slight tendency toward
underestimation. Range, diffusion, decay, and emissions variance were fairly
well identified, while the two variance ratios were poorly identified. The
unidentifiability of the regression variance ratio had little effect on the
regression coefficients however. These three regression coefficients (shown in
Figure 6) were somewhat underestimated, likely due to oversmoothing in the
kriging estimates.

Due to the underestimation issues in the regression coefficients, we recommend
interpreting model output largely in terms of overall trends and predictions
of the source/solution functions. For example, in Figure 5, the model
successfully identifies that the left-most zone has generally higher emissions
than two zones on the right. In cases where accurate estimation of regression
coefficients is a high priority, we recommend specifying non-Gaussian priors
such as the Laplace or Horseshoe priors, and updating the coefficients in the
MCMC step. Ultimately, the model allows us to infer the solution with high
accuracy, and although the source is more uncertain, the inference remains
accurate, with the true source falling well within the range of variation
among the samples.

###  4.2 1-Dimensional Spatio-Temporal Source Results

Figure 7 shows the results of our analysis when we add a temporal dimension to
the 1-D system analyzed in Numerical Case Study 3.6. We simulated the
advection-diffusion equation with a spatio-temporal Matérn source as defined
in equation (15). We simulated 10 sensors equidistant along the $s$-axis, and
simulated 20 measurements each, equidistant in time.

The solution (concentration) estimate is more accurate than the source
(emissions) estimate, because this is the scale at which we are able to take
direct measurements. The source estimate displays a greater degree of
smoothing due to the obscuring effects of diffusion and advection. Because the
material detected in a concentration measurement could have originated from
any point nearby or upstream, a larger number of measurements are required to
provide equivalent resolution on the source function. As more measurements are
recorded, the reconstruction converges to the simulated source. Our method can
be thought of as a filter that treats concentration measurements as an
integral over previous, mostly upstream sources. That is, a measurement of
concentration at time $t$ only provides information about the source process
up to time $t$, and provides much stronger information from upstream locations
(due primarily to advection) than from downstream locations (due to
diffusion). Information about future or downstream source values is primarily
obtained from the correlation structure of the source process. Figure 7 panel
(f) demonstrates this behavior on its margins, where downstream and future
uncertainties are higher than upstream and past uncertainties.

Figure 7: The results of Case Study 1. The source process was spatio-temporal
Matérn as defined in Section 3.4.1 with parameters:
$\tau=2,\alpha=4,\kappa=1$. Flow rate was constant in time, and varied in
space as shown in figure 2. Red points denote sampling locations. A spatial
buffer zone of length 5 on each side is not included in the plot. No temporal
buffer zone was removed.

###  4.3 Air Pollution Results

The results of this analysis can be visualized in Figure 8 and Table 2.
Uncertainty was quantified using the square root of the diagonal elements of
the posterior covariance matrix (i.e., the pointwise posterior standard
deviations). Posterior uncertainty was lower around the main landmass of the
continental US, because this is where the data were collected, and regions of
high estimated emissions and concentration were clustered around some regions
of high population including Texas, California, Florida, and the Pacific
Northwest. Some anomalies are present in the results, notably some rare
negative predictions for PM 2.5 concentrations. These are non-physical, but
occur almost entirely outside the sampling region, likely due to the
artificial oscillatory behavior induced by the FEM discretization of
advection. These could possibly be reduced or removed by using Petrov-Galerkin
elements or a full spatio-temporal model, but we leave these possibilities for
future work.

Parameter | Post. Mean | Post. Std. Dev. | Prior | Prior Quantiles (.025,.975)
---|---|---|---|---
$\rho$ | 456 km | 38 km | $\alpha=4,\beta=0.1$ | (10,80)
$D$ | 8.2km2/hr | 2.8km2/hr | $\alpha=8.5,\beta=1$ | (3.6,14.4)
$r$ | 0.77/hr | 0.42/hr | $\alpha=1.36,\beta=2.94$ | (0.027,1.5)
$\sigma_{f}^{2}$ | 32 $(\mu$g/m${}^{3})^{2}$ | 44 $(\mu$g/m${}^{3})^{2}$ | $\alpha=1.1,\beta=3.9$ | (1,100)
$v$ | 2.1 | 1.7 | $\alpha=1.1,\beta=0.13$ | (0.33,30)
$\sigma_{\epsilon}^{2}$ | 29.6 $(\mu$g/m${}^{3})^{2}$ | 1.7 $(\mu$g/m${}^{3})^{2}$ | NA | NA
$\sigma_{f}$ | 5.0 $\mu$g/m3 | 2.7 $\mu$g/m3 | NA | NA
$\sigma_{\epsilon}$ | 5.4 $\mu$g/m3 | 0.16 $\mu$g/m3 | NA | NA
Table 2: Posterior mean and standard deviations for each parameter of the
advection-diffusion-reaction source model of U.S. air pollution. Also the
priors, and prior quantiles for those parameters that had them. The bottom
three parameters are generated quantities.

The units of emissions in Case Study 3.6.2 are $\mu$g/m3/hour, because
concentration is measured as a density, in units of $\mu$g/m3. Emissions
measured in $\mu$g/hour would be more intuitive, but would require making a
homogeneity assumption to account for the vertical distribution of PM2.5 (we
use this approach to approximate a prior distribution for $\sigma^{2}_{f}$),
or directly modeling vertical distribution with a fully 3-D model. Such a
model is not difficult to mesh using FEniCS but is much larger and more
computationally intensive, and we leave this extension for future work.

Figure 8: The main results of Case Study 3.6.2, parameterized with the NOMADs
climate model (Rutledge et al., 2006), and using EPA PM 2.5 measurements (US
Environmental Protection Agency, 2021). Advection field (a; upper left);
advection speed (b; upper right); PM 2.5 map (c; middle left); PM 2.5 map
posterior standard deviation (d; middle right); PM 2.5 emissions map (e; lower
left); PM 2.5 source map posterior standard deviation (f; lower right). The
buffer region (approximately 4 degrees wide) is set to zero to avoid showing
boundary anomalies. Dots represent locations of data points.

##  5 Conclusion

We presented a method for SPDE-based inference on the source of a linear PDE.
We demonstrated the use of this method for reconstructing the source of an
advection-diffusion-reaction process in a 1-D non-temporal setting, a 1-D
spatio-temporal setting, and in a 2-D non-temporal case study of air pollution
in the United States. We found that the method worked well to estimate
concentrations, exhibiting close to a standard $1/\sqrt{N}$ convergence rate.
We also found that the method worked well to estimate emissions, but with a
reduced error rate that we hypothesized to be due to the destruction of
information by diffusion processes. When applied in a case study of realistic
complexity, we found that reasonable predictions were made, but we
occasionally observed anomalies such as negative concentration. Although these
were limited in spatial extent and had a negligible impact on the overall
analysis, we hypothesize that these artifacts could be reduced or eliminated
in future work by using Petrov-Galerkin elements, or a full spatio-temporal
model. We provide a framework through which our model can be extended to
space-time, and through which efficient likelihood calculations can be
conducted in space-time, but due to the challenges of conducting the
marginalization step in space-time (specifically, maintaining sparsity when
updating the Cholesky factor of $\bm{Q}_{\bm{u}}$ to add in
$\frac{n}{\sigma^{2}_{\epsilon}}\bm{A}^{T}\bm{A}$), leave further spatio-
temporal modeling to future work.

Possible applications of our method include inferring the source of
environmental DNA measurements, non-point source pollution, or stable
isotopes. Other important linear PDEs include Maxwell’s equations and the wave
equation. Using these equations, our method could be applied to MRI
measurements such as in Bao et al. (2002) or earthquake reconstruction such as
in Das & Suhadolc (1996). Our approach uses an SPDE method similar to Lindgren
et al. (2011) to solve the non-point source problem in a Bayesian framework,
providing the opportunity to explicitly incorporate linear regression models
on the source function. This framework also allows us to incorporate
covariance functions for the source, and to calculate uncertainties, in
contrast to many existing inverse problems approaches, which often operate
under either a prediction-only framework (e.g., Ayvaz 2010), or under more
complex frameworks such as the ensemble Kalman filter (e.g., Stroud et al.
2010). By incorporating our methodology and source-inference perspective into
existing hierarchical frameworks (Wikle & Hooten, 2010), process uncertainty
and non-Gaussian observations could also be accommodated, although the use of
a link-function would warp the mechanistic interpretations we seek to preserve
in this study. Our method has a variety of possible applications, and we have
shown that it is both accurate, and implementable using reproducible general-
purpose FEM software (Alnæs et al., 2015) that can easily accommodate
arbitrary changes to the underlying PDE structure.

## References

  * (1)
  * Alnæs et al. (2015) Alnæs, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M. E. & Wells, G. N. (2015), ‘The FEniCS project version 1.5’, Archive of Numerical Software 3(100).
  * Alves et al. (2009) Alves, C., Kress, R. & Serranho, P. (2009), ‘Iterative and range test methods for an inverse source problem for acoustic waves’, Inverse Problems 25(5), 055005.
  * Ayvaz (2010) Ayvaz, M. T. (2010), ‘A linked simulation–optimization model for solving the unknown groundwater pollution source identification problems’, Journal of Contaminant Hydrology 117(1-4), 46–59.
  * Bakka (2018) Bakka, H. (2018), ‘How to solve the stochastic partial differential equation that gives a Matérn random field using the finite element method’, arXiv preprint arXiv:1803.03765 .
  * Bakka et al. (2020) Bakka, H., Krainski, E., Bolin, D., Rue, H. & Lindgren, F. (2020), ‘The diffusion-based extension of the Matérn field to space-time’, arXiv preprint arXiv:2006.04917 .
  * Bao et al. (2002) Bao, G., Ammari, H. & Fleming, J. L. (2002), ‘An inverse source problem for Maxwell’s equations in magnetoencephalography’, SIAM Journal on Applied Mathematics 62(4), 1369–1382.
  * Brenner et al. (2008) Brenner, S. C., Scott, L. R. & Scott, L. R. (2008), The mathematical theory of finite element methods, Vol. 3, Springer.
  * Brooks & Hughes (1982) Brooks, A. N. & Hughes, T. J. (1982), ‘Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations’, Computer Methods in Applied Mechanics and Engineering 32(1-3), 199–259.
  * Byun & Schere (2006) Byun, D. & Schere, K. L. (2006), ‘Review of the governing equations, computational algorithms, and other components of the Models-3 Community Multiscale Air Quality (CMAQ) modeling system’, Applied Mechanics Reviews 59(2), 51–77.
  * Cameletti et al. (2019) Cameletti, M., Gómez-Rubio, V. & Blangiardo, M. (2019), ‘Bayesian modelling for spatially misaligned health and air pollution data through the INLA-SPDE approach’, Spatial Statistics 31, 100353.
  * Cao et al. (2013) Cao, J., Chow, J. C., Lee, F. S., Watson, J. G. et al. (2013), ‘Evolution of PM2.5 measurements and standards in the US and future perspectives for China’, Aerosol and Air Quality Research 13(4), 1197–1211.
  * Chen et al. (2008) Chen, Y., Davis, T. A., Hager, W. W. & Rajamanickam, S. (2008), ‘Algorithm 887: Cholmod, supernodal sparse cholesky factorization and update/downdate’, ACM Transactions on Mathematical Software (TOMS) 35(3), 1–14.
  * Das & Suhadolc (1996) Das, S. & Suhadolc, P. (1996), ‘On the inverse problem for earthquake rupture: The Haskell-type source model’, Journal of Geophysical Research: Solid Earth 101(B3), 5725–5738.
  * Deiner et al. (2017) Deiner, K., Bik, H. M., Mächler, E., Seymour, M., Lacoursière-Roussel, A., Altermatt, F., Creer, S., Bista, I., Lodge, D. M., de Vere, N. et al. (2017), ‘Environmental DNA metabarcoding: Transforming how we survey animal and plant communities’, Molecular Ecology 26(21), 5872–5895.
  * Doi et al. (2017) Doi, H., Inui, R., Akamatsu, Y., Kanno, K., Yamanaka, H., Takahara, T. & Minamoto, T. (2017), ‘Environmental DNA analysis for estimating the abundance and biomass of stream fish’, Freshwater Biology 62(1), 30–39.
  * Doi et al. (2015) Doi, H., Uchii, K., Takahara, T., Matsuhashi, S., Yamanaka, H. & Minamoto, T. (2015), ‘Use of droplet digital pcr for estimation of fish abundance and biomass in environmental DNA surveys’, PloS one 10(3), e0122763.
  * El Badia & Ha-Duong (2002) El Badia, A. & Ha-Duong, T. (2002), ‘On an inverse source problem for the heat equation. application to a pollution detection problem’, Journal of inverse and ill-posed problems 10(6), 585–599.
  * El Badia et al. (2005) El Badia, A., Ha-Duong, T. & Hamdi, A. (2005), ‘Identification of a point source in a linear advection–dispersion–reaction equation: Application to a pollution source problem’, Inverse Problems 21(3), 1121.
  * El Badia & Hamdi (2007) El Badia, A. & Hamdi, A. (2007), ‘Inverse source problem in an advection-dispersion-reaction system: Application to water pollution’, Inverse Problems 23(5), 2103.
  * Engl et al. (1996) Engl, H. W., Hanke, M. & Neubauer, A. (1996), Regularization of Inverse Problems, Vol. 375, Springer Science & Business Media.
  * Grech et al. (2008) Grech, R., Cassar, T., Muscat, J., Camilleri, K. P., Fabri, S. G., Zervakis, M., Xanthopoulos, P., Sakkalis, V. & Vanrumste, B. (2008), ‘Review on solving the inverse problem in EEG source analysis’, Journal of Neuroengineering and Rehabilitation 5(1), 25.
  * Hao et al. (2019) Hao, S., Li, F., Li, Y., Gu, C., Zhang, Q., Qiao, Y., Jiao, L. & Zhu, N. (2019), ‘Stable isotope evidence for identifying the recharge mechanisms of precipitation, surface water, and groundwater in the ebinur lake basin’, Science of the Total Environment 657, 1041–1050.
  * Harrison et al. (2019) Harrison, J. B., Sunday, J. M. & Rogers, S. M. (2019), ‘Predicting the fate of edna in the environment and implications for studying biodiversity’, Proceedings of the Royal Society B 286(1915), 20191409.
  * Lindgren et al. (2011) Lindgren, F., Rue, H. & Lindström, J. (2011), ‘An explicit link between Gaussian fields and Gaussian Markov random fields: The stochastic partial differential equation approach’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73(4), 423–498.
  * Lushi & Stockie (2010) Lushi, E. & Stockie, J. M. (2010), ‘An inverse Gaussian plume approach for estimating atmospheric pollutant emissions from multiple point sources’, Atmospheric Environment 44(8), 1097–1107.
  * Monk et al. (2003) Monk, P. et al. (2003), Finite element methods for Maxwell’s equations, Oxford University Press.
  * Moraga et al. (2017) Moraga, P., Cramb, S. M., Mengersen, K. L. & Pagano, M. (2017), ‘A geostatistical model for combined analysis of point-level and area-level data using INLA and SPDE’, Spatial Statistics 21, 27–41.
  * Panofsky (1969) Panofsky, H. A. (1969), ‘Air pollution meteorology’, American Scientist 57(2), 269–285.
  * Peters et al. (2001) Peters, T. M., Norris, G. A., Vanderpool, R. W., Gemmill, D. B., Wiener, R. W., Murdoch, R. W., Mcelroy, F. F. & Pitchford, M. (2001), ‘Field performance of PM2.5 federal reference method samplers’, Aerosol Science & Technology 34(5), 433–443.
  * Porter & Devaney (1982) Porter, R. P. & Devaney, A. J. (1982), ‘Holography and the inverse source problem’, Journal of the Optical Society of America 72(3), 327–330.
  * Rachel et al. (2008) Rachel, B.-J., Pearson, T. E., Ramos, F. C., Grimes, C. B. & Bruce MacFarlane, R. (2008), ‘Tracking natal origins of salmon using isotopes, otoliths, and landscape geology’, Limnology and Oceanography 53(4), 1633–1642.
  * Rees et al. (2014) Rees, H. C., Maddison, B. C., Middleditch, D. J., Patmore, J. R. & Gough, K. C. (2014), ‘The detection of aquatic animal species using environmental DNA–a review of eDNA as a survey tool in ecology’, Journal of applied ecology 51(5), 1450–1459.
  * Rourke et al. (2022) Rourke, M. L., Fowler, A. M., Hughes, J. M., Broadhurst, M. K., DiBattista, J. D., Fielder, S., Wilkes Walburn, J. & Furlan, E. M. (2022), ‘Environmental DNA (eDNA) as a tool for assessing fish biomass: A review of approaches and future considerations for resource surveys’, Environmental DNA 4(1), 9–33.
  * Rue & Held (2005) Rue, H. & Held, L. (2005), Gaussian Markov Random Fields: Theory and Applications, CRC Press.
  * Rutledge et al. (2006) Rutledge, G. K., Alpert, J. & Ebisuzaki, W. (2006), ‘NOMADS: A climate and weather model archive at the National Oceanic and Atmospheric Administration’, Bulletin of the American Meteorological Society 87(3), 327–342.
  * Sigrist et al. (2012) Sigrist, F., Künsch, H. R. & Stahel, W. A. (2012), ‘A dynamic nonstationary spatio-temporal model for short term prediction of precipitation’, The Annals of Applied Statistics pp. 1452–1477.
  * Sigrist et al. (2015) Sigrist, F., Künsch, H. R. & Stahel, W. A. (2015), ‘Stochastic partial differential equation based modelling of large space–time data sets’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 77(1), 3–33.
  * Song et al. (2018) Song, W., Jia, H., Li, Z. & Tang, D. (2018), ‘Using geographical semi-variogram method to quantify the difference between NO2 and PM2.5 spatial distribution characteristics in urban areas’, Science of the Total Environment 631, 688–694.
  * Stockie (2011) Stockie, J. M. (2011), ‘The mathematics of atmospheric dispersion modeling’, SIAM Review 53(2), 349–372.
  * Stroud et al. (2010) Stroud, J. R., Stein, M. L., Lesht, B. M., Schwab, D. J. & Beletsky, D. (2010), ‘An ensemble Kalman filter and smoother for satellite data assimilation’, Journal of the American Statistical Association 105(491), 978–990.
  * Tucker (2000) Tucker, W. G. (2000), ‘An overview of PM2.5 sources and control strategies’, Fuel Processing Technology 65, 379–392.
  * US Environmental Protection Agency (2021) US Environmental Protection Agency (2021), ‘Air Quality System Data Mart’, Internet Database available at http://www.epa.gov/ttn/airs/aqsdatamart.  Accessed: 2021-02-24.
  * Ver Hoef et al. (2006) Ver Hoef, J. M., Peterson, E. & Theobald, D. (2006), ‘Spatial statistical models that use flow and stream distance’, Environmental and Ecological Statistics 13(4), 449–464.
  * Watson et al. (2000) Watson, J. G., Chow, J. C. & Pace, T. G. (2000), ‘Fugitive dust emissions’, Crops 3(14), 7.
  * Wikle & Hooten (2010) Wikle, C. K. & Hooten, M. B. (2010), ‘A general science-based framework for dynamical spatio-temporal models’, Test 19(3), 417–451.
  * Xepapadeas (2011) Xepapadeas, A. (2011), ‘The economics of non-point-source pollution’, Annual Review of Resource Economics 3(1), 355–373.
  * Yan et al. (2008) Yan, L., Fu, C.-L. & Yang, F.-L. (2008), ‘The method of fundamental solutions for the inverse heat source problem’, Engineering Analysis with Boundary Elements 32(3), 216–222.