# The SPDE approach for spatio-temporal datasets with advection and diffusion

Lucia Clarotto, Denis Allard, Thomas Romary, Nicolas Desassis

► To cite this version:

Lucia Clarotto, Denis Allard, Thomas Romary, Nicolas Desassis. The SPDE approach for spatio-temporal datasets with advection and diffusion. Spatial Statistics, 2024, 62, pp.100847. 10.1016/j.spasta.2024.100847. hal-04663802

HAL Id: hal-04663802

https://hal.inrae.fr/hal-04663802v1

Submitted on 29 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ELSEVIER

Contents lists available at ScienceDirect

Spatial Statistics

journal homepage: www.elsevier.com/locate/spasta

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Check for
updates

# The SPDE approach for spatio-temporal datasets with advection and diffusion

Lucia Clarotto<sup>a,*</sup>, Denis Allard<sup>b,1</sup>, Thomas Romary<sup>c,1</sup>, Nicolas Desassis<sup>c,1</sup>

## ARTICLE INFO

### Keywords:

Spatio-temporal statistics

Stochastic Partial Differential Equations

advection–diffusion

Geostatistics

Solar radiation

## ABSTRACT

In the task of predicting spatio-temporal fields in environmental science using statistical methods, introducing statistical models inspired by the physics of the underlying phenomena that are numerically efficient is of growing interest. Large space-time datasets call for new numerical methods to efficiently process them. The Stochastic Partial Differential Equation (SPDE) approach has proven to be effective for the estimation and the prediction in a spatial context. We present here the advection–diffusion SPDE with first–order derivative in time which defines a large class of nonseparable spatio-temporal models. A Gaussian Markov random field approximation of the solution to the SPDE is built by discretizing the temporal derivative with a finite difference method (implicit Euler) and by solving the spatial SPDE with a finite element method (continuous Galerkin) at each time step. The “Streamline Diffusion” stabilization technique is introduced when the advection term dominates the diffusion. Computationally efficient methods are proposed to estimate the parameters of the SPDE and to predict the spatio-temporal field by kriging, as well as to perform conditional simulations. The approach is applied to a solar radiation dataset. Its advantages and limitations are discussed.

## Contents

<table><tr><td>1.</td><td>Introduction</td><td>2</td></tr><tr><td>2.</td><td>The spatio-temporal advection-diffusion SPDE and its discretization</td><td>3</td></tr><tr><td>2.1.</td><td>Background</td><td>3</td></tr><tr><td>2.2.</td><td>The spatio-temporal advection-diffusion SPDE</td><td>4</td></tr><tr><td>2.3.</td><td>Discretization</td><td>5</td></tr><tr><td>2.4.</td><td>Stabilization of advection-dominated SPDE</td><td>6</td></tr><tr><td>2.5.</td><td>Spatio-temporal Gaussian Markov random field approximation</td><td>6</td></tr><tr><td>3.</td><td>Estimation, prediction and simulation</td><td>7</td></tr><tr><td>3.1.</td><td>Estimation of the parameters</td><td>8</td></tr><tr><td>3.2.</td><td>Prediction by kriging</td><td>9</td></tr><tr><td>3.3.</td><td>Conditional simulations</td><td>10</td></tr><tr><td>3.4.</td><td>Simulation study</td><td>10</td></tr><tr><td>4.</td><td>Application to a solar radiation dataset</td><td>10</td></tr><tr><td>4.1.</td><td>Estimation and prediction</td><td>12</td></tr></table>

<footnote><sup>a</sup> UMR MIA Paris-Saclay, AgroParisTech, INRAE, Université Paris-Saclay, Palaiseau, 91120, France</footnote>

<footnote><sup>b</sup> Biostatistiques et Processus Spatiaux (BioSP), INRAE, Avignon, 84914, France</footnote>

<footnote><sup>c</sup> Centre for geosciences and geoengineering, Mines Paris, PSL University, Fontainebleau, 77300, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: lucia.clarotto@agroparistech.fr (L. Clarotto).</footnote>

<footnote n="1">¹ Contributed equally to this work.</footnote>
<table><tr><td>4.2. Conditional simulations</td><td>14</td></tr><tr><td>5. Conclusion</td><td>14</td></tr><tr><td>Acknowledgments</td><td>16</td></tr><tr><td>Appendix A. Relation between advection-diffusion SPDE and infinite dimensional SDE</td><td>16</td></tr><tr><td>Appendix B. Proof of Proposition 1</td><td>16</td></tr><tr><td>Appendix C. Discretization of spatio-temporal advection-diffusion SPDE</td><td>17</td></tr><tr><td>Appendix D. Proof of Proposition 6</td><td>19</td></tr><tr><td>References</td><td>19</td></tr></table>

# 1. Introduction

Many areas of environmental science seek to predict a space-time variable of interest from observations at scattered points in the space cross time domain of study, e.g., among other possible applications, wind prediction (Lenzi and Genton, 2020; Huang et al., 2022), precipitation forecasting (Sigrist et al., 2011), urban air quality inference (Paciorek et al., 2009). Among modern techniques proposing efficient methods for estimation and prediction in a spatio-temporal framework, there is a distinction between two possible ways of constructing and treating spatio-temporal models (Wikle and Hooten, 2010): either one follows the traditional geostatistical paradigm, using joint space-time covariance functions (see for example Cressie and Huang (1999), Gneiting (2002), Stein (2005), as well as the recent reviews Porcu et al. (2021), Chen et al. (2021)), or one uses dynamical models, including functional time series of surfaces, see for example Wikle and Cressie (1999), Sigrist et al. (2012) and Martínez-Hernández and Genton (2023).

While the theoretical aspects of spatio-temporal geostatistics are well developed (Cressie and Wikle, 2011), their implementation faces difficulties. The geostatistical paradigm is computationally expensive for large spatio-temporal datasets, due to the factorization of dense covariance matrices, whose complexity scales with the cube of the number of observations. This well known problem is often referred to as the “big n problem” (Banerjee et al., 2014). Separable space-time covariance functions have often been used to take advantage of their computational convenience, even when they are not realistic in describing the processes due to the absence of space-time interaction. In most applications, separable models show poorer predictions than nonseparable models, see references above. Recent studies have focused on constructing nonseparable models, which are physically more realistic, albeit computationally more expensive, see Gneiting (2002), Porcu et al. (2006), Salvaña and Genton (2021) and Bourotte et al. (2016), Allard et al. (2022) in a multivariate context. Nonseparable space-time covariance models can be constructed from Fourier transforms of permissible spectral densities, mixtures of separable models, and partial differential equations (PDEs) representing physical laws (Carrizo-Vergara et al., 2022; Lindgren et al., 2022). They can be fully symmetric or asymmetric, stationary or nonstationary, univariate or multivariate, in the Euclidean space or on the sphere. See Porcu et al. (2021) and Chen et al. (2021) for recent comprehensive reviews.

In this paper, we follow the dynamic approach that makes use of physical laws and study models which are defined through Stochastic Partial Differential Equations (SPDEs), where the stochasticity is obtained by adding a random noise as a forcing term. The SPDE approach relies on the representation of a continuously indexed Gaussian Random Field (GRF) as a discretely indexed random process, i.e. a Gaussian Markov Random Field (GMRF, see Rue and Held (2005)). Passing from a GRF to a GMRF, the covariance function and the dense covariance matrix are substituted respectively by a neighborhood structure and a sparse precision matrix. Using GMRFs with sparse precision matrices implies computationally efficient numerical methods, especially for matrix factorization. The link between GRF and GMRFs in the purely spatial case has been pioneered by Lindgren et al. (2011), who proposed to construct a GMRF representation of the spatial Matérn field on a triangulated mesh of the domain through the discretization of a diffusion SPDE with a Finite Element Method (FEM). We refer to Bakka (2022) for a simple explanation of the FEM applied to the spatial SPDE and to Section 2.3 for a detailed generalization to the spatio-temporal SPDE.

In the spatial framework, major mathematical and algorithmic advances in the SPDE approach have been made (Fuglstad et al., 2015; Pereira and Desassis, 2019; Pereira et al., 2022), making it possible to efficiently process very large datasets, even in the presence of nonstationarities and varying local anisotropies. The development of SPDE-based approaches to Gaussian processes has led to several practical solutions, among which we find the R package for approximate Bayesian inference R-INLA (Rue et al., 2009; Lindgren and Rue, 2015) that uses SPDEs to sample from spatial and spatio-temporal models.

When generalizing to the spatio-temporal framework, a direct space-time formulation of the SPDE approach was first suggested in Lindgren et al. (2011), without any precise detail on estimation and prediction. In Cameletti et al. (2011), the SPDE approach was coupled with an AR(1) model in time, leading to a separable space-time model. Nonseparable spatio-temporal models have been elaborated in Särkkä et al. (2013), Krainski et al. (2018) and Lindgren et al. (2023) as a spatio-temporal generalization of the diffusion-Matérn model of Lindgren et al. (2011). In the approaches overviewed above, the space-time processes are symmetrical in the sense that the spatio-temporal covariance does not change when the sign of the space and/or time lag changes. However, atmospheric and geophysical processes are often asymmetric due to transport effects, such as air and water flows. Sigrist et al. (2015) built non-symmetrical and nonseparable space-time Gaussian models as a solution to an advection-diffusion SPDE with computationally efficient algorithms for statistical estimation using fast Fourier transforms and Kalman filters. Sharrock and Kantas (2022) used a similar method, but in an joint online parameter estimation and optimal sensor placement problem. Liu et al. (2020) extended this approach to spatially-varying advection-diffusion and non-zero mean source-sink, leading to a space-time covariance which is nonstationary in space. The applicability of these approaches remains difficult however, especially with scattered data, as it relies on the Fourier transform of the data. Carrizo-Vergara et al. (2022) defined new spatio-temporal models incorporating
the physical processes linked to the studied phenomena (advection, diffusion, etc.), but the estimation of the parameters and the conditioning to the observed data remained unaddressed.

In this work, we propose a new and efficient approach for dealing with spatio-temporal SPDEs that includes both a diffusion and an advection term. In contrast to Sigrist et al. (2015) and Liu et al. (2020), we make use of the sparse formulation of the spatio-temporal field which is the approximate solution of the SPDE obtained by a combination of FEM and finite differences (FD). This sparse formulation allows fast algorithms for parameter estimation and spatio-temporal prediction. We also treat the case of an advection-dominated SPDE, by introducing the streamline diffusion stabilization term in the SPDE (Hughes and Brooks, 1979). To the best of our knowledge, this work is the first statistical FEM/FD implementation of spatio-temporal SPDEs with advection.

The paper is organized as follows: Section 2 first presents background material on the spatio-temporal SPDE approach. The spatio-temporal advection-diffusion model developed in this paper is presented, along with its discretization. Moreover, the stabilization of advection-dominated SPDEs is introduced. Section 3 explores fast and scalable estimation methods, kriging formula for prediction and conditional simulations. Section 4 presents an application of the proposed spatio-temporal SPDE approach to a solar radiation dataset. Section 5 discusses the advantages and the limitations of the approach and opens the way to further works.

## 2. The spatio-temporal advection-diffusion SPDE and its discretization

### 2.1. Background

In the SPDE representation, GRFs on $\mathbb{R}^d$ are viewed as solutions to specific stochastic partial differential equations (Whittle, 1954, 1963). In particular, Gaussian Whittle-Matérn fields, analyzed in details in Lindgren et al. (2011) and reviewed in Lindgren et al. (2022), are solutions to

$$ (\kappa^2 - \Delta)^{\alpha/2} X(\cdot) = \tau W(\cdot), \qquad (1) $$

with $\alpha > d/2$ and $\tau > 0$. $\Delta = \sum_{i=1}^{d} \frac{\partial^2}{\partial s_i^2}$ is the Laplacian operator and $W(\cdot)$ is a standard spatial Gaussian white noise, whose definition is briefly recalled.

A white noise $W(\cdot)$ is as a Generalized Random Field (GeRF) that associates to any function $\phi \in L^2(\mathbb{R}^d)$ a random variable $W(\phi) \in \mathbb{R}$, that satisfies

$$ E[W(\phi)] = 0, \quad \forall \phi \in L^2(\mathbb{R}^d) $$

and

$$ \mathrm{Cov}[W(\phi_1), W(\phi_2)] = \int_{\mathbb{R}^d} \phi_1(s)\phi_2(s) \, ds, \quad \forall \phi_1, \phi_2 \in L^2(\mathbb{R}^d). \qquad (2) $$

If, moreover, for any $m \ge 1$ and any linearly independent $\phi_1, \dots, \phi_m \in L^2(\mathbb{R}^d)$, the random vector $[W(\phi_1), \dots, W(\phi_m)]^T$ is a Gaussian vector, then $W$ is called Gaussian white noise.

If $\{Z_i\}_{i \in \mathbb{N}}$ is a sequence of independent, standard Gaussian variables, then the function $W$ defined over $L^2(\mathbb{R}^d)$ by

$$ W(\phi) = \sum_{j \in \mathbb{N}} Z_j \int_{\mathbb{R}^d} \phi e_j \, ds, \quad \forall \phi \in L^2(\mathbb{R}^d), $$

where $\{e_j\}_{j \in \mathbb{N}}$ denotes an orthonormal basis of $L^2(\mathbb{R}^d)$, is a Gaussian white noise on $\mathbb{R}^d$.

In principle, GeRFs have only meaning when applied to test functions in some particular functional space, and not necessarily when evaluated in points of the space, but, for an easier reading, we will allow ourselves to write $W(s)$ and $X(s)$.

The covariance function of the Gaussian Whittle-Matérn field solution to Eq. (1) is the well known Matérn covariance function

$$ \mathrm{Cov}(\mathbf{h}) = \sigma^2 C_{\nu}^{M}(\kappa \|\mathbf{h}\|) = \frac{\sigma^2}{2^{\nu-1} \Gamma(\nu)} (\kappa \|\mathbf{h}\|)^{\nu} \mathcal{K}_{\nu}(\kappa \|\mathbf{h}\|), \qquad (3) $$

with smoothness parameter $\nu = \alpha - d/2 > 0$, scale parameter $\kappa$ and variance $\sigma^2 = \tau^2(4\pi)^{-d/2}\Gamma(\nu)\Gamma(\nu + d/2)^{-1}\kappa^{-2\nu}$. $\mathcal{K}_{\nu}$ is the modified 2nd order Bessel function and $\mathbf{h} = \mathbf{s} - \mathbf{s}'$ is the spatial lag between two locations $\mathbf{s}$ and $\mathbf{s}'$ in $\mathbb{R}^d$. In particular, when $\nu = 1/2$, we get the exponential covariance function and when $\nu \to +\infty$, after proper renormalization, (3) tends to the Gaussian covariance function (Genton, 2001).

In Lindgren et al. (2011), the smoothness parameter $\nu$ considered in the Matérn covariance function corresponds to integer values of $\alpha$. When non-integer values of $\alpha$ are introduced in the modeling, the SPDE is said to be fractional. Recent reviews of results and applications of the fractional SPDE approach are available in Bolin and Kirchner (2020), Roques et al. (2022), Bolin et al. (2024), but this case will not be treated further in this work.

When generalizing to spatio-temporal processes $X(t,s)$, we consider the framework proposed in Carrizo-Vergara et al. (2022) for extending the SPDE approach to a wide class of linear spatio-temporal SPDEs. Let us denote $\xi \in \mathbb{R}^d$ a spatial frequency and $\omega \in \mathbb{R}$ a temporal frequency. The space-time white noise with unit variance, denoted $W(t,s)$, is characterized by its spectral measure $d\mu_W(\omega, \xi) = (2\pi)^{-(d+1)}d\xi d\omega$. New spatio-temporal models were obtained from known PDEs describing physical processes, such as diffusion, advection, and oscillations with stochastic forcing terms. In particular, Carrizo-Vergara et al. (2022) provided sufficient conditions to the existence and uniqueness of stationary solutions to

$$ \left[ \frac{\partial \beta}{\partial t^{\beta}} + \mathcal{L}_g \right] X(t, s) = W(t, s), \qquad (4) $$
with $\beta > 0$. In (4), the spatial operator $\mathcal{L}_g$ is defined using the spatial Fourier transform on $\mathbb{R}^d$, denoted $\mathcal{F}_S$,

$$ \mathcal{L}_g(\cdot) = \mathcal{F}_S^{-1}(g \mathcal{F}_S(\cdot)), $$

where $g: \mathbb{R}^d \to \mathbb{C}$ is a sufficiently regular and Hermitian-symmetric function called the symbol function of the operator $\mathcal{L}_g$. The temporal operator $\frac{\partial^\beta}{\partial t^\beta}$ is

$$ \frac{\partial^\beta}{\partial t^\beta}(\cdot) = \mathcal{F}_T^{-1}((i\omega)^\beta \mathcal{F}_T(\cdot)), $$

where $\mathcal{F}_T$ is the temporal Fourier transform on $\mathbb{R}$ and where we have used the symbol function over $\mathbb{R}$.

$$ \omega \mapsto (i\omega)^\beta = |\omega|^\beta e^{i\operatorname{sgn}(\omega)\beta\pi/2}. $$

The spatio-temporal symbol function of the operator involved in (4) is thus

$$ (\omega, \xi) \mapsto (i\omega)^\beta + g(\xi) = |\omega|^\beta \cos\left(\frac{\beta\pi}{2}\right) + g_R(\xi) + i \left(\operatorname{sgn}(\omega)|\omega|^\beta \sin\left(\frac{\beta\pi}{2}\right) + g_I(\xi)\right), $$

where $g_R$ and $g_I$ are the real and imaginary part of the spatial symbol function $g(\xi)$. If $|g_R|$ is inferiorly bounded by the inverse of a strictly positive polynomial and $g_R \cos(\frac{\beta\pi}{2}) \ge 0$, Theorem 1 and Proposition 3 in Carrizo-Vergara et al. (2022) state that (4) admits a unique stationary solution for every arbitrary $g_I$ function.

## 2.2. The spatio-temporal advection–diffusion SPDE

The advection–diffusion equation is a Partial Differential Equation (PDE) that describes physical phenomena where particles, energy, or other physical quantities evolve inside a physical system due to two processes: diffusion and advection. Advection represents the mass transport due to the average velocity of all particles, and diffusion represents the mass transport due to the instantaneously varying velocity of individual particles. In this paper, we study the advection–diffusion SPDE on the domain $[0, T] \times \mathbb{R}^d$ that writes

$$ \left[ \frac{\partial}{\partial t} + \frac{1}{c}(\kappa^2 - \nabla \cdot \mathbf{H} \nabla)^\alpha + \frac{1}{c}\gamma \cdot \nabla \right] X(t,s) = \frac{\tau}{\sqrt{c}} Z(t,s), \quad (5) $$

where

* the operator $\nabla \cdot \mathbf{H} \nabla$ is a diffusion term that can incorporate anisotropy in the matrix $\mathbf{H}$. When the field is isotropic, i.e. when $\mathbf{H} = \lambda \mathbf{I}$, this term reduces to the Laplacian operator $\lambda \Delta$;

* the operator $\gamma \cdot \nabla$ models the advection, $\gamma \in \mathbb{R}^d$ being a velocity vector;

* $\alpha \ge 0$ relates to the smoothness of $X(t, \cdot)$, $\kappa^2 > 0$ accounts for damping and $c$ is a positive time-scale parameter;

* $\tau \ge 0$ is a standard deviation factor and $Z$ is a stochastic forcing term, detailed below.

This equation was mentioned in Lindgren et al. (2011), Carrizo-Vergara et al. (2022) and Lindgren et al. (2023), and was analyzed using spectral approaches in Sigrist et al. (2015) and Liu et al. (2020). The term $Z(t,s)$ is assumed to be of the form

$$ Z(t,s) = W_T(t) \otimes Z_S(s), \quad (6) $$

i.e., a space–time separable stochastic (generalized) random function given as the tensor product of a temporal Gaussian white noise $W_T$ and a spatial noise $Z_S$. $Z_S$ is often chosen to be a spatial Gaussian white noise, denoted $W_S$ in this case. To ensure a sufficient regularity for $Z$, $Z_S$ can alternatively be a colored noise, such as for example the solution to the spatial Whittle–Matérn SPDE (Lindgren et al., 2011)

$$ (\kappa^2 - \nabla \cdot \mathbf{H} \nabla)^{\alpha/2} Z_S(s) = W_S(s), \quad (7) $$

where $W_S$ is a Gaussian white noise. Note that the parameters $\kappa^2$ and $\mathbf{H}$ in the noise term have been set identical to those in the diffusion term in the left-hand-side of (5) to ensure that the spatial marginalization of the process is a Matérn field, as detailed below. A relation can be found between the SPDE notation of our paper and the more classical notation of infinite dimensional SDEs of Da Prato and Zabczyk (1992). This relation is explained in Appendix A.

When $\alpha > 0$, $X(t,s)$ is a stationary nonseparable spatio-temporal field with covariance function $\operatorname{Cov}(u, h)$, with $(u, h) \in \mathbb{R} \times \mathbb{R}^d$. The advection–diffusion Eq. (5) is a particular first–order evolution model as in Eq. (4) with $\beta = 1$. Its spatial symbol function,

$$ g(\xi) = \frac{1}{c} [(\kappa^2 + \xi^\top \mathbf{H} \xi)^\alpha + i\gamma^\top \xi], $$

verifies the sufficient condition for existence and uniqueness of a stationary solution recalled at the end of Section 2.1. We define the spatial trace of $X$ as the spatial random field $X(t, \cdot)$ at any $t \in [0, T]$. Carrizo-Vergara et al. (2022) showed that the advection term does not affect the spatial trace of the solution. For some specific values of the parameters, the spatial trace of the solution to (5) is a Matérn field, as detailed in Proposition 1. In the following $|\mathbf{H}|$ denotes the determinant of the square matrix $\mathbf{H}$.
**Proposition 1.** Let $Z(t,s)$ be a spatio-temporal noise colored in space with $Z_S(s)$ satisfying (7), and let $\alpha_{tot} = \alpha + \alpha_S$. If $\alpha_{tot} > d/2$, the spatial trace of the stationary solution $X(t,s)$ of the SPDE (5) is the Gaussian Matérn field with covariance

$$ \mathrm{Cov}(0, \mathbf{h}) = \frac{\tau^2 \Gamma(\alpha_{tot} - d/2)}{2\Gamma(\alpha_{tot})(4\pi)^{d/2}\kappa^{2(\alpha_{tot}-d/2)}|\mathbf{H}|^{1/2}} C_{\alpha_{tot}-d/2}^{M} \left( \kappa \|\mathbf{H}^{-1/2}\mathbf{h}\| \right), \quad (8) $$

where $\mathbf{h} = \mathbf{s} - \mathbf{s}'$ is the spatial lag and $C_{\alpha_{tot}-d/2}^{M}(\cdot)$ is the unit variance and scale Matérn covariance function defined in (3) with smoothness parameter equal to $\nu = \alpha_{tot} - d/2$.

**Proposition 1** is adapted from Proposition 1 in Lindgren et al. (2023). A proof is reported in Appendix B. The model reduces to a separable one in a particular case stated in the corollary below.

**Corollary 2.** Let the coefficients of the SPDE (5) be such that $\alpha = 0$ and $\gamma = 0$; the spatial operator applied to the spatio-temporal field $X(t,s)$ is then the constant value $c^{-1}$. Let $Z(t,s)$ be a spatio-temporal noise colored in space, with $Z_S(s)$ satisfying (7). If $\alpha_S > d/2$, the stationary solution of the SPDE is a separable spatio-temporal field with covariance

$$ \mathrm{Cov}(u, \mathbf{h}) = \frac{\tau^2 \Gamma(\alpha_S - d/2)}{2\Gamma(\alpha_S)(4\pi)^{d/2}\kappa^{2(\alpha_S - d/2)}|\mathbf{H}|^{1/2}} C_{\alpha_S - d/2}^{M} \left( \kappa \|\mathbf{H}^{-1/2}\mathbf{h}\| \right) \exp \frac{u}{c}, $$

with smoothness parameter equal to $\nu = \alpha_S - d/2$.

## 2.3. Discretization

The advection-diffusion SPDE in (5) is discretized in time and space, using finite differences and a finite element method, respectively (from now on, this type of discretization will be noted as FEM/FD discretization). The temporal domain $[0, T]$ is discretized in $(N_T + 1)$ regular time steps of length $dt = T/N_T$, and we note $t_k = kdt$ for $k \in \{0, \dots, N_T\}$. Since implicit solvers are usually less sensitive to numerical instability than explicit solvers, the implicit Euler scheme is chosen for the temporal discretization. This choice implies stability, hence convergence towards the stationary solution. We denote $X^{(k)} = X(t_k, \cdot)$ the temporal approximation of the spatial trace of the SPDE (5) at time $t_k$. The FEM method for the spatial discretization is the continuous Galerkin method with Neumann Boundary Conditions as detailed in Lindgren et al. (2011).

The solution in two dimensions is now detailed. The solution in three dimensions involve geometrical technicalities, but is otherwise very similar. Let $\Omega \subset \mathbb{R}^2$ be a compact and connected domain of $\mathbb{R}^2$. $\Omega$ is meshed using a triangulation $\mathcal{T}$ with $N_S$ vertices $\{s_1, \dots, s_{N_S}\} \subset \Omega$. Let $h := \max_{\mathrm{Tr} \in \mathcal{T}} h_{\mathrm{Tr}}$, where $h_{\mathrm{Tr}}$ is the length of the longest side of the triangle $\mathrm{Tr} \in \mathcal{T}$. A first-order finite element representation $X_h$ of the solution to the spatial SPDE is a linear combination $X_h = \sum_{i=1}^{N_S} x_i \psi_i$ of piecewise linear basis functions $\{\psi_i\}_{i=1}^{N_S}$, each $\psi_i$ being equal to 1 at the vertex $s_i$ and 0 at all the other vertices. The weights $\{x_i\}_{i=1}^{N_S}$ define uniquely the values of the field at the vertices, while the values in the interior of the triangles are determined by linear interpolation. The continuous Galerkin solution is then obtained by finding the weights that fulfill the weak formulation of Eq. (5) for test functions belonging to the space $\mathcal{V}$ spanned by $\{\psi_i\}_{i=1}^{N_S}$.

**Proposition 3.** Let $X(t,s)$ be the spatio-temporal process solution to Eq. (5) with $\alpha \in \{0, 1\}$ and spatio-temporal white noise, i.e. $Z(t,s) = W(t,s) = W_T(t) \otimes W_S(s)$. Let $\mathcal{T}$ be a triangulation of $\Omega$ and $\{\psi_i\}_{i=1}^{N_S}$ be the piecewise linear basis functions defined over $\mathcal{T}$. Let us define the mass matrix $\mathbf{M} = [M_{ij}]_{i,j=1}^{N_S}$, the stiffness matrix $\mathbf{G} = [G_{ij}]_{i,j=1}^{N_S}$, the advection matrix $\mathbf{B} = [B_{ij}]_{i,j=1}^{N_S}$ and the matrix $\mathbf{K} = [K_{ij}]_{i,j=1}^{N_S}$ as follows:

$$ M_{ij} = \int_{\Omega} \psi_i(s) \psi_j(s) \mathrm{d}s, $$

$$ G_{ij} = \int_{\Omega} \mathbf{H} \nabla \psi_i(s) \cdot \nabla \psi_j(s) \mathrm{d}s, $$

$$ B_{ij} = \int_{\Omega} \gamma \cdot \nabla \psi_i(s) \psi_j(s) \mathrm{d}s, $$

$$ K_{ij} = (\kappa^2 M_{ij} + G_{ij})^a. $$

Then, at each time step, the continuous Galerkin finite element solution vector $\mathbf{x}^{(k+1)} = \{x_i^{(k+1)}\}_{i=1}^{N_S}$, for $k \in \{0, \dots, N_T\}$, satisfies

$$ \left( \mathbf{M} + \frac{dt}{c} (\mathbf{K} + \mathbf{B}) \right) \mathbf{x}^{(k+1)} = \mathbf{M} \mathbf{x}^{(k)} + \frac{\tau \sqrt{dt}}{\sqrt{c}} \mathbf{M}^{1/2} \mathbf{z}^{(k+1)}, \quad (9) $$

where $\mathbf{z}^{(k+1)} \sim \mathcal{N}(0, \mathbf{I}_{N_S})$, $\mathbf{M}^{1/2}$ is any matrix such that $\mathbf{M}^{1/2}\mathbf{M}^{1/2} = \mathbf{M}$ and $dt = T/N_T$. When the noise on the right-hand side is colored in space, i.e. $Z(t,s) = W_T(t) \otimes Z_S(s)$, the discretization reads

$$ \left( \mathbf{M} + \frac{dt}{c} (\mathbf{K} + \mathbf{B}) \right) \mathbf{x}^{(k+1)} = \mathbf{M} \mathbf{x}^{(k)} + \frac{\tau \sqrt{dt}}{\sqrt{c}} \mathbf{M} \mathbf{L}_S^{\top} \mathbf{z}^{(k+1)}, $$

where $\mathbf{L}_S$ is the Cholesky decomposition of $\mathbf{Q}_S^{-1}$, the covariance matrix of the discretized solution $\mathbf{Z}_S$ of the spatial SPDE (7), obtained with the continuous Galerkin FEM (Lindgren et al., 2011).
**Proof.** The proof is available in Appendix C.

**Remark 1.** The elements of the matrices M, G and B are non-zero only for pairs of basis functions which share common triangles. This implies that the matrix $(M + \frac{dt}{c}(K + B))$ is sparse and that Eq. (9) can be solved by Cholesky decomposition in an efficient way.

## 2.4. Stabilization of advection-dominated SPDE

When the advection term is too strong with respect to the diffusion term, advection-domination occurs. In the framework outlined above, when $\alpha = 1$, the non-symmetric matrix $[M + \frac{dt}{c}(K + B)]$ becomes ill-conditioned, which induces oscillations and unstable solutions for the continuous Galerkin approximation. Specifically, the advection-domination occurs when the Péclet number $Pe^h = \frac{\|\gamma\|_h}{2\lambda} > 1$, where $\lambda$ is the coefficient of the isotropic Laplacian operator (see for example Mekuria and Rao (2016) or Quarteroni (2008, Chapter 5)).

One possible solution is to decrease the mesh size $h$, i.e., to refine the triangulation, until the advection no longer dominates on the element-level, with $Pe^h < 1$. However, in many cases this is not a feasible solution because it would increase the number of vertices beyond computation limits. Another solution, adopted here, is to introduce a stabilization term. Many stabilization approaches are possible, some being more accurate than others (Quarteroni, 2008, Chapter 5). In our case, we opt for the streamline diffusion stabilization approach (Hughes and Brooks, 1979), considered as a good trade-off between accuracy and computational complexity. Essentially, the SD approach consists of stabilizing the advection by introducing an artificial diffusion term along the advection direction. The following proposition presents the stabilized solution to (5).

**Proposition 4.** Assume the same hypotheses as in Proposition 3 with $\alpha = 1$. The solution to Eq. (5) in presence of streamline diffusion stabilization is

$$ (M + \frac{dt}{c}(K + B + S)) x^{(k+1)} = M x^{(k)} + \frac{\bar{\tau} \sqrt{dt}}{\sqrt{c}} M^{1/2} z^{(k+1)}, \quad (10) $$

where $S = [S_{ij}]_{i,j=1}^{NS}$ is the matrix of the streamline diffusion stabilization operator S, such that

$$ S_{ij} = S(\psi_i, \psi_j) = h \|\gamma\|^{-1} \int_{\Omega} (\gamma \cdot \nabla \psi_i)(\gamma \cdot \nabla \psi_j) \, d\mathbf{s}, $$

and $\bar{\tau} = \tau (|\mathbf{H} + h\|\gamma\|^{-1}\gamma\gamma^\top|)^{-1/4} (|\mathbf{H}|)^{1/4}$. When the noise on the right-hand side of Eq. (5) is colored in space, the discretization becomes

$$ (M + \frac{dt}{c}(K + B + S)) x^{(k+1)} = M x^{(k)} + \frac{\bar{\tau} \sqrt{dt}}{\sqrt{c}} M L_S^\top z^{(k+1)}, $$

where $L_S$ is as in Proposition 3.

The proof of the discretized equation follows the same reasoning as that of Proposition 3 with the addition of the matrix S. The SD approach can be seen as a perturbation of the original SPDE (Bank et al., 1990). Indeed, by making the classical hypothesis of Neumann boundary condition on $\Omega$ and by using the Green's first identity, we get

$$ \int_{\Omega} (\gamma \cdot \nabla x)(\gamma \cdot \nabla v) \, d\mathbf{s} = - \int_{\Omega} \nabla \cdot (\gamma \gamma^\top) \nabla x v \, d\mathbf{s}. $$

As a consequence, the original SPDE (5) can be rewritten with an additional diffusion term as

$$ \left[ \frac{\partial}{\partial t} + \frac{1}{c} \left[ \kappa^2 - \nabla \cdot (\mathbf{H} + h \|\gamma\|^{-1} \gamma \gamma^\top) \nabla + \gamma \cdot \nabla \right] \right] X(t, s) = \frac{\tau}{\sqrt{c}} Z(t, s). \quad (11) $$

The term $(h\|\gamma\|^{-1}\gamma\gamma^\top)$ acts as an anisotropic “diffusion” matrix that is added to the anisotropy (or identity) matrix $\mathbf{H}$ of the original diffusion. This extra diffusion stabilizes the advection directed along the direction $\gamma$. By following the proof of Proposition 1, we find that the marginal variance of the spatial field $X(t, \cdot)$ of Eq. (11) is equal to

$$ \sigma^2 = \frac{\tau^2 \Gamma(\alpha_{\text{tot}} - d/2)}{\Gamma(\alpha_{\text{tot}}) 2(4\pi)^{d/2} \kappa^{2(\alpha_{\text{tot}} - d/2)} |\mathbf{H} + h \|\gamma\|^{-1} \gamma \gamma^\top|^{1/2}}. $$

For the variance to be equal to the variance in Proposition 1, $\tau$ must be replaced by $\bar{\tau} = \tau (|\mathbf{H} + h\|\gamma\|^{-1}\gamma\gamma^\top|)^{1/4} (|\mathbf{H}|)^{-1/4}$.

## 2.5. Spatio-temporal Gaussian Markov random field approximation

**Proposition 5.** In presence of an advection-dominated flow and a spatio-temporal white noise on the right-hand side of Eq. (5), the discretized vector $x^{(k+1)}$ on the mesh $\mathcal{T}$ at each time step is the solution of the following equation:

$$ \begin{gathered} x^{(0)} \sim \mathcal{N}(0, \Sigma), \\ x^{(k+1)} = \mathbf{D}x^{(k)} + \mathbf{E}z^{(k+1)}, \quad (12) \end{gathered} $$
where

$$
\begin{align}
\mathbf{D} &= \left( \mathbf{M} + \frac{dt}{c} (\mathbf{K} + \mathbf{B} + \mathbf{S}) \right)^{-1} \mathbf{M}, \\
\mathbf{E} &= \frac{\bar{\tau} \sqrt{dt}}{\sqrt{c}} \left( \mathbf{M} + \frac{dt}{c} (\mathbf{K} + \mathbf{B} + \mathbf{S}) \right)^{-1} \mathbf{M}^{1/2}, \tag{13}
\end{align}
$$

and $\mathbf{x}^{(k+1)} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}_{N_S})$ is independent of $\mathbf{x}^{(0)}, \dots, \mathbf{x}^{(k+1)}$. In presence of a spatio-temporal noise colored in space on the right-hand side of Eq. (5), the matrix $\mathbf{E}$ reads

$$
\mathbf{E} = \frac{\bar{\tau} \sqrt{dt}}{\sqrt{c}} \left( \mathbf{M} + \frac{dt}{c} (\mathbf{K} + \mathbf{B} + \mathbf{S}) \right)^{-1} \mathbf{M} \mathbf{L}_S^{\top},
$$

where $\mathbf{L}_S$ is defined in Proposition 3.

**Proof.** Starting from Eq. (10), which represents the numerical scheme for the advection–diffusion spatio-temporal SPDE with stabilization, it is straightforward to obtain (12).

When the SPDE is not advection-dominated, which implies that no stabilization term is needed, Eq. (13) is replaced by the similar equation where the matrix S is deleted and $\bar{\tau}$ is replaced by $\tau$.

The covariance matrix $\Sigma$ of the approximate spatial trace $\mathbf{x}^{(0)}$ at the first time step, can be taken to be equal to any admissible positive definite matrix. The closer $\Sigma$ is to the covariance $C_S$ of $X(t, \cdot)$, the faster the stationary solution is obtained. When the hypotheses of Proposition 1 are satisfied, an efficient option is to choose $\Sigma$ as the Matérn covariance of Eq. (8).

To obtain fast inference and prediction computations, the precision matrix of the spatio-temporal discretized solution $\mathbf{x}_{0:N_T} = [\mathbf{x}^{(0)}, \dots, \mathbf{x}^{(N_T)}]^{\top}$ must be sparse. For this reason, $\mathbf{M}$ is replaced by the diagonal matrix $\tilde{\mathbf{M}} = [\tilde{M}_{ij}]_{i,j=1}^{N_S}$, where $\tilde{M}_{ii} = \langle \psi_i, 1 \rangle$ and $\tilde{M}_{ij} = 0$ if $i \neq j$ (Lindgren et al., 2011). This technique is called mass lumping and is common practice in FEM (Quarteroni, 2008, Chapter 5). From now on, we always use the diagonal matrix $\tilde{\mathbf{M}}$, but for ease of reading, it will still be denoted $\mathbf{M}$.

**Proposition 6.** Let $\mathbf{x}_{0:N_T} = [\mathbf{x}^{(0)}, \dots, \mathbf{x}^{(N_T)}]^{\top}$ be the vector containing all spatial solutions until time step $N_T$ of Eq. (12). The global precision matrix $\mathbf{Q}$ of the vector $\mathbf{x}_{0:N_T}$ of size $(N_S(N_T + 1), N_S(N_T + 1))$ reads

$$
\mathbf{Q} = \begin{pmatrix}
\Sigma^{-1} + \mathbf{D}^{\top}\mathbf{F}^{-1}\mathbf{D} & -\mathbf{D}^{\top}\mathbf{F}^{-1} & 0 & \cdots & 0 \\
-\mathbf{F}^{-1}\mathbf{D} & \mathbf{F}^{-1} + \mathbf{D}^{\top}\mathbf{F}^{-1}\mathbf{D} & -\mathbf{D}^{\top}\mathbf{F}^{-1} & \cdots & \mathbf{F}^{-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & -\mathbf{F}^{-1}\mathbf{D} & \mathbf{F}^{-1} + \mathbf{D}^{\top}\mathbf{F}^{-1}\mathbf{D} & -\mathbf{D}^{\top}\mathbf{F}^{-1} & 0 \\
0 & \cdots & 0 & -\mathbf{D}^{\top}\mathbf{F}^{-1} & \mathbf{F}^{-1}
\end{pmatrix}, \quad (14)
$$

where $\mathbf{F} = \mathbf{E}\mathbf{E}^{\top}$.

**Proof.** The proof is available in Appendix D.

**Remark 2.** Matrix (14) has a tridiagonal structure in time and is sparse in each of its spatial blocks of size $(N_S, N_S)$ located on the three diagonals. The sparsity of the precision matrix relies on the choice of $\tilde{\mathbf{M}}$ and on the value of $\alpha$ (the higher $\alpha$ the less sparse the spatial blocks). Conversely, the precision of finite element method discretization is influenced by the mesh density (the smaller $h$, the more precise the solution); this factor plays a role in defining the size of the precision matrix, but not its sparsity pattern. The sparsity pattern will be useful for the following sections concerning estimation and prediction.

3. Estimation, prediction and simulation

This section presents an efficient implementation for parameter estimation and spatio-temporal prediction within the spatio-
temporal SPDE framework described in Section 2. We consider the advection–diffusion SPDE (5) with $d = 2$, $\alpha = 1$, $\mathbf{H} = \mathbf{I}$ (isotropic
diffusion) and colored noise in space with $\alpha_S = 2$. Similar computations can be generalized to other values of $\alpha_S$ such that $\alpha_S/2$ is
integer or to anisotropic diffusion.

The spatio-temporal domain $\Omega \times [0, T]$ is discretized in space with a triangulation $\mathcal{T}$ with $N_S$ nodes and discretized in time
by means of $(N_T + 1)$ regular time steps. This space-time discretization is denoted $\mathcal{T}' = \mathcal{T} \times \{0, T/N_T, \dots, T\}$. At each time step
$k \in \{0, \dots, N_T\}$ there are $n^{(k)}$ observations scattered in the spatial domain $\Omega$. There is thus a total of $n = \sum_{k=0}^{N_T} n^{(k)}$ spatio-temporal
data collected in the vector $\mathbf{y}_{0:N_T} = [\mathbf{y}^{(0)}]^{\top}, \dots, [\mathbf{y}^{(N_T)}]^{\top}]^{\top}$.

We consider a statistical model with fixed and random effects. The fixed effect is a linear trend on a set of covariates and the random effect is modeled as the FEM/FD discretization of the random field solution to the SPDE (5) (as described in Section 2.3), with the addition of random noise:

$$
\mathbf{y}_{0:N_T} = \mathbf{\Lambda}\mathbf{b} + \mathbf{\Lambda}\mathbf{x}_{0:N_T} + \sigma_0\varepsilon, \quad (15)
$$

where **b** is the vector of *q* fixed effects and **Λ** is a (n, q) matrix of covariates with [Λ]_{jk} = λ_k(t_j, s_j), j = 1..., n and k = 1, ..., q.
The matrix **A** is the (n, N_S(N_T + 1)) projection matrix between the data and the points in T', and ε is a standard Gaussian
random vector with independent components. When the observation locations do not change during the time window, A^T A is a $(N_S(N_T + 1), N_S(N_T + 1))$ block-diagonal matrix with all $(N_S, N_S)$ equal blocks.

The discretization with FEM in space and FD in time is justified by a few considerations: compared to a Fourier approximation approach (such as the ones considered in Sigrist et al. (2015), Liu et al. (2020), Sharrock and Kantas (2022)), it is better suited for scattered data, it can be easily adapted to spatially and/or temporally varying parameters in the SPDE, leading to nonstationary extensions of the method, and it can be generalized to spatio-temporal fields on Riemannian surfaces with only a few changes in the approach (Pereira, 2023).

### 3.1. Estimation of the parameters

The parameters of the SPDE are estimated using maximum likelihood. We collect the parameters of the SPDE in the vector $\theta^T = [\kappa, \gamma_x, \gamma_y, c, \tau]$, while all the parameters of the statistical model are collected in $\psi^T = [\theta^T, b^T, \sigma_0]$. Following (15), $y_{0:N_T}$ is a Gaussian vector with expectation $\Lambda b$ and covariance matrix

$$ \Sigma_{y_0:N_T} = \Lambda Q^{-1}(\theta)\Lambda^T + \sigma_0^2 I_n, $$

where $Q(\theta)$ is a precision matrix of size $(N_S(N_T + 1), N_S(N_T + 1))$ depending on the parameters $\theta$. For ease of notation, we use $Q$ instead of $Q(\theta)$. The log-likelihood is equal to

$$ \mathcal{L}(\psi) = -\frac{n}{2}\log(2\pi) - \frac{1}{2}\log|\Sigma_{y_0:N_T}| - \frac{1}{2}(y_{0:N_T} - \Lambda b)^T \Sigma_{y_0:N_T}^{-1} (y_{0:N_T} - \Lambda b). \quad (16) $$

We use the Broyden, Fletcher, Goldfarb, and Shanno optimization algorithm (Nocedal and Wright, 2006), that makes use of the second-order derivative of the objective function. The gradients of the log-likelihood function (16) with respect to the different parameters included in $\psi$ are approximately computed with finite differences.

We now propose and detail a computationally efficient formulation for both the terms of the log-likelihood (16), namely $-\frac{1}{2}\log|\Sigma_{y_0:N_T}|$ and $-\frac{1}{2}(y_{0:N_T} - \Lambda b)^T \Sigma_{y_0:N_T}^{-1} (y_{0:N_T} - \Lambda b)$. We first consider the log-determinant. The quadratic form is addressed next.

**Proposition 7.** In the framework outlined above, we have

$$ \log|\Sigma_{y_0:N_T}| = n \log \sigma_0^2 - \log|Q| + \log|Q_A|, \quad (17) $$

where $Q_A = Q + \sigma_0^{-2} \Lambda^T A$.

**Proof.** To compute $\log|\Sigma_{y_0:N_T}|$, let us consider the augmented matrix

$$ \Sigma_c = \begin{pmatrix} Q^{-1} & Q^{-1} \Lambda^T \\ \Lambda Q^{-1} & \Sigma_{y_0:N_T} \end{pmatrix}. \quad (18) $$

Hence,

$$ Q_c = \Sigma_c^{-1} = \begin{pmatrix} Q + \sigma_0^{-2} \Lambda^T A & -\sigma_0^{-2} \Lambda^T \\ -\sigma_0^{-2} A & \sigma_0^{-2} I_n \end{pmatrix}. \quad (19) $$

Using block formulas, we have

$$ \log|\Sigma_c| = -\log|Q_c| = -\log|Q| + n \log \sigma_0^2, $$

and

$$ \log|\Sigma_c| = \log|\Sigma_{y_0:N_T}| + \log|Q^{-1} - Q^{-1}\Lambda^T \Sigma_{y_0:N_T}^{-1} \Lambda Q^{-1}| \\ = \log|\Sigma_{y_0:N_T}| - \log|Q + \sigma_0^{-2} \Lambda^T A|, $$

where the last equality is a consequence of the Woodbury identity. This leads to the result.

**Proposition 8.** The term $\log|Q|$ in Eq. (17) can be computed with the computationally efficient formula

$$ \log|Q| = \log|\Sigma^{-1}| + (N_T - 1) \log|F^{-1}|, \quad (20) $$

with

$$ F^{-1} = \frac{c}{\varepsilon^2 dt} (M + \frac{dt}{c}(K + B + S))^T M^{-1} Q_S M^{-1} (M + \frac{dt}{c}(K + B + S)), $$

where $Q_S$ is the precision matrix of the discretized spatial noise $Z_S$ defined in Proposition 3.
**Proof.** Following Powell (2011), let $N_N = [N_{ij}]_{i,j=1}^N$ be an $(nN, nN)$ matrix, which is partitioned into $N^2$ blocks $N_{ij}$, each of size $(n, n)$. Then the determinant of $N_N$ is

$$|N_N| = \prod_{k=1}^{N} |\alpha_{kk}^{(N-k)}|,$$

where the $\alpha^{(k)}$ are defined by

$$\alpha_{ij}^{(0)} = N_{ij} \\ \alpha_{ij}^{(k+1)} = \alpha_{ij}^{(k)} - \alpha_{i,N-k}^{(k)} \alpha_{N-k,N-k}^{(k)}^{-1} \alpha_{N-k,j}^{(k)}, \quad k \ge 1.$$

Q is a block-matrix organized as $N_N$. Hence, the formula for $|Q|$ is

$$|Q| = |\Sigma^{-1}||F^{-1}|^{(N_T+1)-1}. \qquad (21)$$

Applying the logarithm, we obtain Eq. (20).

Note that $|F^{-1}|$ is now the determinant of a $(N_S, N_S)$ sparse, symmetric and positive definite matrix. The log-determinant can be computed through its Cholesky decomposition as

$$\log|F^{-1}| = 2\log\left(\prod_{i=1}^{n} F_{\text{chol},ii}^{-1}\right) = 2 \sum_{i=1}^{N_S} \log(F_{\text{chol},ii}^{-1}),$$

since $F_{\text{chol}}^{-1}$ is a triangular matrix, whose determinant is the product of the diagonal elements.

The term $\log|Q_A| = \log|Q + \sigma_0^{-2}A^T A|$ requires a detailed analysis. The term $\sigma_0^{-2}A^T A$ is an $(N_S(N_T + 1), N_S(N_T + 1))$ diagonal block matrix, whose $(N_S, N_S)$ blocks are sparse. The computation of $\log|Q_A|$ is not as straightforward as in the case of $\log|Q|$, because there is no way of reducing the computation to purely spatial matrices. Depending on the size $N_S(N_T + 1)$, we can either apply a Cholesky decomposition of the $(N_S(N_T + 1), N_S(N_T + 1))$ matrix $Q_A$ or the matrix-free approach proposed in Pereira et al. (2022), here sketched. The logarithm function is first approximated by a Chebyshev polynomial $P(\cdot)$ (Chebyshev, 1853), then the Hutchinson's estimator (Hutchinson, 1990) is used to obtain a stochastic estimate of tr[$P(Q_A)$]. The method is detailed in Algorithm 5 in Pereira et al. (2022).

Now, concerning the quadratic term of the log-likelihood (16), we can work with the more convenient expression obtained thanks to the Woodbury formula

$$\Sigma_{y_0:N_T}^{-1} = \sigma_0^{-2} I_n - \sigma_0^{-4} A Q_A^{-1} A^T.$$

Hence

$$(y_{0:N_T} - \mathbf{A} \mathbf{b})^T \Sigma_{y_0:N_T}^{-1} (y_{0:N_T} - \mathbf{A} \mathbf{b}) = \sigma_0^{-2} (y_{0:N_T} - \mathbf{A} \mathbf{b})^T \mathbf{I}_n (y_{0:N_T} - \mathbf{A} \mathbf{b}) - \sigma_0^{-4} (y_{0:N_T} - \mathbf{A} \mathbf{b})^T \mathbf{A} \mathbf{Q}_A^{-1} \mathbf{A}^T (y_{0:N_T} - \mathbf{A} \mathbf{b}). \quad (22)$$

The second term of Eq. (22) can be computed either by Cholesky decomposition or using the conjugate gradient (CG) method. This latter method solves $Q_A v = w$ with respect to $v$ and computes $v_{\text{sol}} = w^T v$, with $w = A^T (y_{0:N_T} - \mathbf{A} \mathbf{b})$. In this case, it is useful to find a good preconditioner for the matrix $Q_A$ to ensure fast convergence of the CG method. We found that a temporal block Gauss-Seidel preconditioner (Young, 1971, Chapter 3) was a good choice in this case. A detailed presentation of the CG method is available in Pereira et al. (2022).

## 3.2. Prediction by kriging

Under a Gaussian assumption, optimal prediction is the conditional expectation, also known in the geostatistics literature as kriging. We detail here two prediction settings: space-time interpolation and temporal extrapolation.

In the space-time interpolation setting, the spatio-temporal vector $x_{0:N_T}$ is predicted on the entire spatial mesh during the time window $[0, T]$, i.e. on $T'$, using the data $y_{0:N_T}$ defined in Eq. (15). The kriging predictor is directly read from Eq. (19):

$$\mathbf{x}_{0:N_T}^* = \mathrm{E}(\mathbf{x}_{0:N_T} | \mathbf{y}_{0:N_T}) = \sigma_0^{-2} \mathbf{Q}_A^{-1} \mathbf{A}^T (\mathbf{y}_{0:N_T} - \mathbf{A} \mathbf{b}). \quad (23)$$

The CG method is used to avoid inverting $Q_A$ in Eq. (23), as detailed at the end of Section 3.1. The conditional variance, also called kriging variance, is

$$\mathrm{Var}(\mathbf{x}_{0:N_T} | \mathbf{y}_{0:N_T}) = \mathbf{Q}_A^{-1}.$$

The computation of the diagonal of an inverse matrix is not straightforward when only the Cholesky decomposition of the matrix is available. Among the existing methods there is the Takahashi recursive algorithm described in Takahashi et al. (1973) and Erisman and Tinney (1975). Another way of computing the kriging variance is through conditional simulations, as detailed in Section 3.3.

In the temporal extrapolation setting, the vector $x^{(N_T+1)}$ is predicted at time step $(N_T + 1)$ on $T$ using all the data available until time $T$, i.e. from $y_{0:N_T}$. Following Eq. (12), we have

$$\mathbf{x}^{(N_T+1)} = \mathbf{D} \mathbf{x}^{(N_T)} + \mathbf{E} \mathbf{z}^{(N_T+1)}, \quad (24)$$
<caption>**Table 1**</caption>

Mean (and standard deviation) of ML estimates $\hat{\theta}^T = [\hat{\kappa}, \hat{\gamma}_1, \hat{\gamma}_2, \hat{\epsilon}, \hat{\tau}]$ over 50 simulations for two different subsets of advection-diffusion model parameters and two different estimation approaches (Cholesky decomposition and matrix-free approach).

<table><thead><tr><th>Method</th><th>κ</th><th>γ<sub>1</sub></th><th>γ<sub>2</sub></th><th>c</th><th>τ</th><th>Average time (s)</th></tr></thead><tbody><tr><td>correct</td><td>0.5</td><td>2</td><td>2</td><td>1</td><td>1</td><td></td></tr><tr><td>Cholesky</td><td>0.610 (0.047)</td><td>2.354 (0.515)</td><td>2.325 (0.421)</td><td>1.037 (0.218)</td><td>1.072 (0.040)</td><td>194</td></tr><tr><td>Matrix-free</td><td>0.483 (0.029)</td><td>1.904 (0.178)</td><td>1.906 (0.147)</td><td>1.027 (0.046)</td><td>0.998 (0.024)</td><td>960</td></tr><tr><td>correct</td><td>0.7</td><td>1</td><td>-1</td><td>2</td><td>0.5</td><td></td></tr><tr><td>Cholesky</td><td>0.695 (0.067)</td><td>1.056 (0.659)</td><td>-1.134 (0.631)</td><td>2.090 (0.352)</td><td>0.503 (0.022)</td><td>172</td></tr><tr><td>Matrix-free</td><td>0.669 (0.052)</td><td>0.967 (0.123)</td><td>-1.164 (0.142)</td><td>1.954 (0.118)</td><td>0.485 (0.013)</td><td>863</td></tr></tbody></table>

where $\mathbf{z}^{(N_T+1)}$ is a standardized Gaussian vector, and $\mathbf{D}$ and $\mathbf{E}$ are defined in Proposition 5. The kriging predictor $\mathbf{x}^{*(N_T+1)}$ is

$$ \mathbf{x}^{*(N_T+1)} = \mathbb{E}(\mathbf{x}^{(N_T+1)} | \mathbf{y}_{0:N_T}) = \mathbf{D} \mathbb{E}(\mathbf{x}^{(N_T)} | \mathbf{y}_{0:N_T}) = \mathbf{D} \mathbf{x}^{*(N_T)}, \quad (25) $$

where $\mathbf{x}^{*(N_T)}$ is extracted from $\mathbf{x}_{0:N_T}^*$. The same procedure can be iterated to predict $\mathbf{x}$ at further time steps.

### 3.3. Conditional simulations

To perform a conditional simulation, we use the conditional kriging paradigm presented below. This approach relies on the fact that kriging predictors and kriging residuals are uncorrelated (independent under Gaussian assumption, see Chilès and Delfiner (1999, Chapter 7)). First, a non-conditional simulation $\mathbf{x}_{0:N_T}^{(NC)}$ is performed on the spatio-temporal grid $\mathcal{T}'$. From this simulation, kriging residuals

$$ \mathbf{r}_{0:N_T} = \mathbb{E} \left( \mathbf{x}_{0:N_T} | \mathbf{A} \mathbf{x}_{0:N_T}^{(NC)} \right) - \mathbf{x}_{0:N_T}^{(NC)} $$

are computed over the entire spatio-temporal grid $\mathcal{T}'$. The conditional expectation is computed using the method presented in the previous section. In a second step, these independently generated residuals are added to the usual kriging of the data to get the conditional simulation

$$ \mathbf{x}_{0:N_T}^{(C)} = \mathbf{x}_{0:N_T}^{*} + \mathbf{r}_{0:N_T}. $$

Conditional simulations at further time steps are obtained by iteratively computing $\mathbf{x}_{N_T+k}^{(C)}$ using the propagation Eq. (24) with $k \ge 1$. Multiple independent realizations of conditional simulations can then be used to compute estimates of conditional variances or other quantities, such as probability maps of threshold exceedance.

### 3.4. Simulation study

We report here some results regarding the estimation of the parameters $\theta^T = [\kappa, \gamma_1, \gamma_2, c, \tau]$ for 50 independent realizations of a spatio-temporal model simulated with the SPDE (5). We set $\mathbf{H} = \mathbf{I}$, $\alpha = 1$ and $\alpha_S = 2$. The spatial domain is the $[0, 30]^2$ square with a grid triangulation of $N_S = 900$ spatial points. The time window is $[1, 10]$ with unit time step and $N_T + 1 = 10$. The $n_S = 100$ observations are randomly located into the spatial domain and their position do not change during the $N_T$ time steps (hence $n = 1000$). Since the sizes of both the dataset and the spatio-temporal mesh are reasonable, we report the estimations computed with both the Cholesky decomposition approach and the matrix-free approach.

As initial values, we use estimated values obtained from the variograms of the spatial and temporal traces of the process. Specifically, the initial value for $\kappa$ is the estimated scale parameter of a Matérn covariance function with smoothness parameter $\nu = \alpha + \alpha_S - 1 = 2$ considering independent temporal repetitions, the initial value for $c$ is deduced from the estimated parameter of $n_S$ independent repetitions of AR(1) processes of length $(N_T+1)$ and $\tau^2$ is computed from Eq. (8) with $\sigma^2$ being the empirical variance computed on the data. Finally, the initial value for $\gamma$ is the null vector. The parameters for the matrix-free approach are set to the following: the order of the Chebyshev polynomial to approximate the logarithm is set to 30 and the number of terms in the sum of the Hutchinson's estimator is set to 10. The results are reported in Table 1. They show that all parameters are accurately estimated with both approaches. In almost all cases, the true value of the parameter is within the mean $\pm 2$ standard deviations interval. We remark how the matrix-free approach takes more time to estimate the parameters. This is due to the iterative computations, that increase the computational time. However, we know that the benefit of the matrix-free approach is the possibility of applying it to much larger spatio-temporal meshes, where the Cholesky decomposition cannot be applied at all.

## 4. Application to a solar radiation dataset

With the constantly increasing installation of photovoltaic (PV) power and its volatility due to weather conditions, characterizing short-term variability of generated solar power (from the minute resolution to the 15-minute resolution) is important for the integration of PV systems into the electrical grid, for balancing supply and demand (Kreuwel et al., 2020). Fluctuations in solar production can have a significant impact on grid stability, and accurate prediction allows for planning necessary adjustments, such as modulating the production of other energy sources, to avoid service interruptions and optimize resource utilization.
<figure><img src="image_5.png" /><figcaption>Fig. 1. Left: spatial domain of study (red square). Right: zoom on the spatial domain along with measuring stations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 2. GHI G and Clear Sky Index K_c for 4 different stations on the 28th of May 2013.</figcaption></figure>

The approach detailed in the previous sections is now applied to a solar radiation dataset for which experts agree on the presence of advection due to Western prevailing winds transporting clouds from one side of the domain to the other. The HOPE campaign (Macke et al., 2017) recorded Global Horizontal Irradiance (GHI) (also called SSI, Surface Solar Irradiance) over a 10 × 16 km² region in West Germany near the city of Jülich from April 2 to July 2, 2013. The sensors were located at 99 stations located as pictured in Fig. 1 and GHI was recorded every 15 s. A detailed description of the campaign can be found in Macke et al. (2017).

The dataset was cleaned for outlying values and non-operating sensors, and the temporal resolution was reduced from 15 s to 1 min. Fig. 2 (left panel) shows GHI as a function of time (in minute, during a full day — the 28th of May 2013) at 4 different stations. These stations, represented in color in Fig. 1, are located at the border of the domain, far from each other. The GHI starts close to 0, increases after sunrise, peaks at midday and tends to 0 at sunset. The maximal theoretical amount of irradiance reaching the sensor follows an ideal concave curve. The divergence between the measured irradiance and the optimal curve can be slight or important, depending on the presence of clouds. One can see on this example that the evolution among the 4 stations is similar, with variations accounting for spatio-temporal variations of the clouds.

A first preprocessing was made in order to stationarize the phenomenon. Oumbe et al. (2014) showed that the solar irradiance at ground level, GHI (denoted G for short from now on), computed by a radiative transfer model can be approximated by the product of the irradiance under clear atmosphere (called Clear Sky GHI, or G_c) and a modification factor due to cloud properties and ground albedo only (Clear Sky Index, or K_c, Beyer et al. (1996)):

$$G \approx G_c K_c. \qquad (26)$$

The error made using this approximation depends mostly on the solar zenith angle, the ground albedo and the cloud optical depth. In most cases, the maximum errors (95th percentile) on global and direct surface irradiances are less than 15 Wm⁻² and less than 2 to 5% in relative value, as recommended by the World Meteorological Organization for high-quality measurements of the solar irradiance (Oumbe et al., 2014). Practically, it means that a model for fast calculation of surface solar irradiance may be separated into two distinct and independent models: a deterministic model for G_c, under clear-sky conditions, as computed according to Gschwind et al. (2019), considered as known in this study, and a model for K_c, which accounts for cloud influence on the downwelling radiation and is expected to change in time and space. K_c is modeled as a random spatio-temporal process and will be the subject of our analysis. Fig. 2 (right panel) shows the variable K_c corresponding to the variable G shown on the left panel. In general, K_c lies between 0 and 1, but in rare occasions, values above 1 can be observed. This phenomenon is called overshooting (Schade et al., 2007) and is due to light reflection by surrounding clouds.
<figure><img src="image_7.png" /><figcaption>Fig. 3. Left: Histogram of $K_c$ over 20 time steps. Right: Time series of $K_c$ for all the stations over 20 time steps, along with the mean value for each time step (in red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<table><caption>**Table 2**
Estimated parameters and log-likelihood for 6 different models from all data on a 20-minute window.</caption><thead><tr><th></th><th>Model</th><th>α</th><th>α<sub>S</sub></th><th>log-likelihood</th><th>κ̂</th><th>γ̂<sub>x</sub></th><th>γ̂<sub>y</sub></th><th>κ̂</th><th>τ̂</th><th>σ̂<sub>0</sub></th><th>μ̂</th></tr></thead><tbody><tr><td>(1)</td><td>adv-diff</td><td>0</td><td>2</td><td>2587</td><td>1.477</td><td>9.642</td><td>-5.382</td><td>11.659</td><td>2.254</td><td>0.052</td><td>0.570</td></tr><tr><td>(2)</td><td>adv-diff</td><td>1</td><td>0</td><td>2577</td><td>0.237</td><td>4.718</td><td>-0.928</td><td>9.315</td><td>0.458</td><td>0.045</td><td>0.598</td></tr><tr><td>(3)</td><td>adv-diff</td><td>1</td><td>2</td><td>2579</td><td>1.299</td><td>17.325</td><td>-8.442</td><td>41.017</td><td>3.072</td><td>0.058</td><td>0.574</td></tr><tr><td>(4)</td><td>diff</td><td>0</td><td>2</td><td>2507</td><td>1.240</td><td>0</td><td>0</td><td>12.558</td><td>1.081</td><td>0.059</td><td>0.569</td></tr><tr><td>(5)</td><td>diff</td><td>1</td><td>0</td><td>2545</td><td>0.246</td><td>0</td><td>0</td><td>6.594</td><td>0.436</td><td>0.047</td><td>0.577</td></tr><tr><td>(6)</td><td>diff</td><td>1</td><td>2</td><td>2512</td><td>0.940</td><td>0</td><td>0</td><td>34.607</td><td>1.248</td><td>0.064</td><td>0.580</td></tr></tbody></table>

A time window of 20 min around 4 p.m. on May 28, 2013 is extracted with observations every minute at the 73 stations with well recorded values. The histogram and time series of the data are shown in Fig. 3. Parameters are estimated on this 20-minute window using the method described in Section 3.1. The spatio-temporal grid contains $N_T + 1 = 20$ one-minute time steps, from $t = 1$ to $t = 20$ and $N_S = 900$ regular spatial mesh points.

## 4.1. Estimation and prediction

Six different models are fitted to the data and used for prediction: 3 models with advection (called “adv-diff”) and 3 models without advection (called “diff”) obtained by setting $γ = 0$. Both groups contain the three following sub-models: (i) a model with diffusion included only in the stochastic forcing term, with a Matérn spatial trace with $ν = 1$; (ii) a nonseparable model that does not have a Matérn spatial trace, but a generalized covariance function instead (the spatial covariance function exists only if $α + α_S > d/2$, where $d = 2$ in this case, see Proposition 1); in this case, the SPDE generates a process which only has meaning as a random measure, and cannot strictly be interpreted at individual locations; (iii) a nonseparable model with a Matérn spatial trace with $ν = 2$. In the general model of Eq. (5) they correspond respectively to $(α, α_S) = (0, 2)$, $(1, 0)$, $(1, 2)$. The parameters of the statistical model of Eq. (15) are estimated for each model separately. The fixed effects term does not include any covariates and is denoted $μ$. The results are reported in Table 2.

The log-likelihoods of the models that include advection are within a range of variations of 10 log-likelihood units and are between 34 to 80 units larger than those with diffusion only. As a point of comparison, if all spatio-temporal dependencies were ignored, the BIC penalization for the advection parameters would be $2\ln(1460) \approx 14.5$. These results indicate strong evidence in favor of models with advection, but no significant differences among them. The parameters vary substantially from one model to the other, but it must be remembered that, when considered independently, their physical interpretation is model dependent. Some combinations are interpretable however. For example, following Proposition 1, the overall variance is equal to $(8\pi)^{-1}\tau^2\kappa^{-2}$ (or $(8\pi)^{-1}\tau^2\kappa^{-2}|\mathbf{I} + h\|\gamma\|^{-1}\gamma\gamma^T|^{-1/2}$ in the stabilized case) when $α_{tot} = α + α_S = 2$ and it is equal to $(16\pi)^{-1}\tau^2\kappa^{-4}$ (or $(16\pi)^{-1}\tau^2\kappa^{-4}|\mathbf{I} + h\|\gamma\|^{-1}\gamma\gamma^T|^{-1/2}$) when $α_{tot} = 3$. Accordingly, the estimated standard deviations for models (1), (3), (4) and (6) are equal to 0.160, 0.119, 0.174 and 0.199 respectively, with the experimental standard deviation being equal to 0.184. For the same models, the practical ranges computed as $\sqrt{8ν/κ}$ (Lindgren et al., 2011) are equal to 1.915, 3.079, 2.281 and 4.255 respectively. Notice that among pairs of models that differ by the presence or absence of advection, the estimated range is larger for those without advection in an attempt to account for the larger correlation distance due to transport.

With the objective of improving the predictions at the 15-minute resolution needed for the electrical grid management, we then perform prediction with two different validation settings containing 80% of conditioning data and 20% of validation data. In the first case (called “Uniform”) the validation locations are uniformly randomly selected. In the second case (called “South-East”) the validation locations are located downwind (i.e. South-East) with respect to the estimated advection direction. See Fig. 4 for a representation of the validation settings.

Recall that a time window containing 20 unit time-steps, from $t = 1$ to $t = 20$, has been selected. For each validation setting and considering each time the end of the time window from $T = 11$ to $T = 20$, three prediction configurations using conditioning data from time $(T - 9)$ to $T$ are computed and compared to the real values, allowing us to compute a root mean square error
<figure><img src="image_8.png" /><figcaption>Fig. 4. Validation settings: Uniform (left) and South-East (right).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 5. Predictions of $K_c$ at (T+1), (T+2) and (T+3) with model (3) ("adv-diff" with $\alpha = 1$ and $\alpha_s = 2$). The black contoured dots are the conditioning locations and the white contoured dots are the validation locations for the Uniform setting. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 6. Averaged RMSE for ST prediction configuration for 6 different models and 2 validation settings: Uniform (left) and South-East (right).</figcaption></figure>

(RMSE) validation score. First, the kriging is performed spatially only (hereafter referred to as “S” kriging). Second, a temporal extrapolation is computed at the conditioning locations at time horizons (T + 1), (T + 2) and (T + 3) (“T1”, “T2”, “T3” kriging). Third, the spatio-temporal prediction is computed at the validation locations at time horizons (T + 1), (T + 2) and (T + 3) (“ST1”, “ST2”, “ST3” kriging). We thus have a total of 6 models × 2 validation settings × 3 prediction configurations. RMSEs are averaged over the 10 repetitions. Results are reported in Table 3 and those of the ST prediction configuration are also shown in Fig. 6.

For all tested validation settings and prediction configurations, the models with advection show better RMSE scores than models without advection. This result is a confirmation of the results already observed on log-likelihoods. Models with advection have similar prediction scores in the prediction configurations S and T, model (1) having slightly better performances in the configuration S. In the *T* and ST configurations, models (2) and (3) have in general quite similar RMSEs, except in the South-East setting with ST configuration where model (2) is clearly the best model. In this case, prediction is made in a space–time domain lying downstream with respect to the advection. It is thus expected that the model best representing the underlying physics should lead to the best prediction performances.

An example of prediction maps on $\mathcal{T}$ at time horizons (T + 1), (T + 2) and (T + 3) is reported in Fig. 5, along with the observed values (black contoured dots). The white contoured dots are the locations used for validation in the Uniform setting.
<caption>**Table 3**</caption>

Averaged RMSE computed at 10 successive time steps for 6 different models, 2 validation settings (Uniform and South-East) and 3 prediction configurations (S, T and ST); see text for details. In each case, the best score among the models is in bold font.

<table><thead><tr><th rowspan="2"></th><th rowspan="2">Model</th><th rowspan="2">α</th><th rowspan="2">α<sub>S</sub></th><th colspan="3">Uniform</th></tr><tr><th>S (min,max)</th><th>T1 (min,max)</th><th>T2 (min,max)</th><th>T3 (min,max)</th></tr></thead><tbody><tr><td>(1)</td><td>adv-diff</td><td>0</td><td>2</td><td></td><td></td><td></td></tr><tr><td>(2)</td><td>adv-diff</td><td>1</td><td>0</td><td></td><td></td><td></td></tr><tr><td>(3)</td><td>adv-diff</td><td>1</td><td>2</td><td></td><td></td><td></td></tr><tr><td>(4)</td><td>diff</td><td>0</td><td>2</td><td></td><td></td><td></td></tr><tr><td>(5)</td><td>diff</td><td>1</td><td>0</td><td></td><td></td><td></td></tr><tr><td>(6)</td><td>diff</td><td>1</td><td>2</td><td></td><td></td><td></td></tr><tr><td></td><td>Model</td><td>α</td><td>α<sub>S</sub></td><td>T1 (min,max)</td><td>T2 (min,max)</td><td>T3 (min,max)</td></tr><tr><td>(1)</td><td>adv-diff</td><td>0</td><td>2</td><td>0.095 (0.067,0.120)</td><td>0.146 (0.111,0.186)</td><td>0.181 (0.131,0.236)</td></tr><tr><td>(2)</td><td>adv-diff</td><td>1</td><td>0</td><td><strong>0.071 (0.046,0.090)</strong></td><td><strong>0.093 (0.054,0.124)</strong></td><td><strong>0.102 (0.060,0.143)</strong></td></tr><tr><td>(3)</td><td>adv-diff</td><td>1</td><td>2</td><td>0.072 (0.046,0.093)</td><td>0.095 (0.054,0.127)</td><td>0.104 (0.055,0.144)</td></tr><tr><td>(4)</td><td>diff</td><td>0</td><td>2</td><td>0.094 (0.058,0.123)</td><td>0.137 (0.091,0.181)</td><td>0.166 (0.102,0.231)</td></tr><tr><td>(5)</td><td>diff</td><td>1</td><td>0</td><td>0.079 (0.058,0.098)</td><td>0.108 (0.082,0.135)</td><td>0.124 (0.099,0.158)</td></tr><tr><td>(6)</td><td>diff</td><td>1</td><td>2</td><td>0.083 (0.054,0.108)</td><td>0.110 (0.077,0.149)</td><td>0.125 (0.085,0.180)</td></tr><tr><td></td><td>Model</td><td>α</td><td>α<sub>S</sub></td><td>ST1 (min,max)</td><td>ST2 (min,max)</td><td>ST3 (min,max)</td></tr><tr><td>(1)</td><td>adv-diff</td><td>0</td><td>2</td><td>0.105 (0.067,0.144)</td><td>0.147 (0.106,0.193)</td><td>0.179 (0.140,0.231)</td></tr><tr><td>(2)</td><td>adv-diff</td><td>1</td><td>0</td><td>0.091 (0.052,0.131)</td><td>0.103 (0.062,0.161)</td><td>0.110 (0.071,0.165)</td></tr><tr><td>(3)</td><td>adv-diff</td><td>1</td><td>2</td><td><strong>0.085 (0.050,0.127)</strong></td><td><strong>0.094 (0.058,0.142)</strong></td><td><strong>0.100 (0.061,0.157)</strong></td></tr><tr><td>(4)</td><td>diff</td><td>0</td><td>2</td><td>0.123 (0.072,0.186)</td><td>0.150 (0.081,0.237)</td><td>0.170 (0.116,0.257)</td></tr><tr><td>(5)</td><td>diff</td><td>1</td><td>0</td><td>0.104 (0.070,0.153)</td><td>0.122 (0.095,0.181)</td><td>0.131 (0.095,0.187)</td></tr><tr><td>(6)</td><td>diff</td><td>1</td><td>2</td><td>0.108 (0.073,0.153)</td><td>0.126 (0.085,0.188)</td><td>0.134 (0.082,0.199)</td></tr><tr><th colspan="7" style="text-align: center;">South-East</th></tr><tr><td></td><td>Model</td><td>α</td><td>α<sub>S</sub></td><td>S (min,max)</td><td>T1 (min,max)</td><td>T2 (min,max)</td></tr><tr><td>(1)</td><td>adv-diff</td><td>0</td><td>2</td><td></td><td></td><td></td></tr><tr><td>(2)</td><td>adv-diff</td><td>1</td><td>0</td><td></td><td></td><td></td></tr><tr><td>(3)</td><td>adv-diff</td><td>1</td><td>2</td><td></td><td></td><td></td></tr><tr><td>(4)</td><td>diff</td><td>0</td><td>2</td><td></td><td></td><td></td></tr><tr><td>(5)</td><td>diff</td><td>1</td><td>0</td><td></td><td></td><td></td></tr><tr><td>(6)</td><td>diff</td><td>1</td><td>2</td><td></td><td></td><td></td></tr><tr><td></td><td>Model</td><td>α</td><td>α<sub>S</sub></td><td>T1 (min,max)</td><td>T2 (min,max)</td><td>T3 (min,max)</td></tr><tr><td>(1)</td><td>adv-diff</td><td>0</td><td>2</td><td>0.095 (0.065,0.122)</td><td>0.142 (0.106,0.185)</td><td>0.172 (0.121,0.228)</td></tr><tr><td>(2)</td><td>adv-diff</td><td>1</td><td>0</td><td>0.074 (0.045,0.099)</td><td>0.097 (0.065,0.128)</td><td>0.109 (0.069,0.148)</td></tr><tr><td>(3)</td><td>adv-diff</td><td>1</td><td>2</td><td><strong>0.074 (0.049,0.097)</strong></td><td><strong>0.096 (0.062,0.122)</strong></td><td><strong>0.106 (0.061,0.139)</strong></td></tr><tr><td>(4)</td><td>diff</td><td>0</td><td>2</td><td>0.090 (0.063,0.116)</td><td>0.128 (0.097,0.167)</td><td>0.154 (0.113,0.209)</td></tr><tr><td>(5)</td><td>diff</td><td>1</td><td>0</td><td>0.081 (0.057,0.105)</td><td>0.111 (0.090,0.147)</td><td>0.128 (0.100,0.169)</td></tr><tr><td>(6)</td><td>diff</td><td>1</td><td>2</td><td>0.084 (0.056,0.109)</td><td>0.109 (0.086,0.152)</td><td>0.123 (0.091,0.176)</td></tr><tr><td></td><td>Model</td><td>α</td><td>α<sub>S</sub></td><td>ST1 (min,max)</td><td>ST2 (min,max)</td><td>ST3 (min,max)</td></tr><tr><td>(1)</td><td>adv-diff</td><td>0</td><td>2</td><td>0.102 (0.081,0.121)</td><td>0.158 (0.130,0.180)</td><td>0.210 (0.160,0.241)</td></tr><tr><td>(2)</td><td>adv-diff</td><td>1</td><td>0</td><td><strong>0.079 (0.044,0.119)</strong></td><td><strong>0.075 (0.039,0.127)</strong></td><td><strong>0.070 (0.038,0.132)</strong></td></tr><tr><td>(3)</td><td>adv-diff</td><td>1</td><td>2</td><td>0.090 (0.052,0.121)</td><td>0.099 (0.056,0.149)</td><td>0.109 (0.057,0.172)</td></tr><tr><td>(4)</td><td>diff</td><td>0</td><td>2</td><td>0.128 (0.092,0.157)</td><td>0.165 (0.112,0.195)</td><td>0.199 (0.116,0.236)</td></tr><tr><td>(5)</td><td>diff</td><td>1</td><td>0</td><td>0.107 (0.050,0.201)</td><td>0.100 (0.042,0.223)</td><td>0.094 (0.0380,0.217)</td></tr><tr><td>(6)</td><td>diff</td><td>1</td><td>2</td><td>0.114 (0.060,0.204)</td><td>0.111 (0.058,0.229)</td><td>0.108 (0.059,0.223)</td></tr></tbody></table>

## 4.2. Conditional simulations

Fig. 7 shows 100 conditional simulations of $K_c$ computed at time $T = 11$ and horizons $(T + 1), (T + 2), \dots, (T + 6)$ with the advection-diffusion model (3). Two validation stations have been selected: one in the North-West part of the domain (the orange star in the left panel of Fig. 4) and one in the South-East part of the domain (the green star in the right panel of Fig. 4). Given that there is an advection from NW to SE, it is therefore expected that the advection-diffusion model should be able to transport the information in that direction. The mean of the 100 simulations and the envelopes corresponding to twice the pointwise standard deviation have also been represented, along with the true values. As expected, most of the conditional simulations lie within the envelopes in both cases and at all time horizons. The remarkable result is that the variance of the conditional simulations at the green station is smaller than that at the orange one at every time step, especially when the time horizon increases. This is due to the advection term in model (3), able to propagate information from North-West to South-East.

## 5. Conclusion

The spatio-temporal SPDE approach based on advection-diffusion equations proposed in this work combines elements of physics, numerical analysis and statistics. It can be seen as a first step towards physics informed geostatistics, which introduces
Acknowledgments

We are grateful to the O.I.E. center of Mines Paris, especially to Yves-Marie Saint-Drenan, Philippe Blanc and Hadrien Verbois, for providing the data and for inspiring discussions about renewable resources evaluation. We thank the Mines Paris/INRAE chair “Geolearning” for its constant support. We are thankful to the STSDS department of King Abdullah University of Science and Technology (KAUST), and especially to Professor Marc Genton, for the insightful work carried out during a stay in Saudi Arabia.

Appendix A. Relation between advection–diffusion SPDE and infinite dimensional SDE

Let us consider the advection–diffusion SPDE (5), here recalled:

$$
\left[ \frac{\partial}{\partial t} + \frac{1}{c} (\kappa^2 - \nabla \cdot \mathbf{H} \nabla)^\alpha + \frac{1}{c} \gamma \cdot \nabla \right] X(t, \mathbf{s}) = \frac{\tau}{\sqrt{c}} Z(t, \mathbf{s}), \qquad (A.1)
$$

$$
Z(t, s) = W_T(t) \otimes Z_S(s). \tag{A.2}
$$

The spatio-temporal forcing as defined in Eq. (A.2) is a GeRF on functions of $L^2([0, T]) \times L^2(\mathbb{R}^d)$. Then, for any $(\phi_T, \phi_S), (\psi_T, \psi_S) \in L^2([0, T]) \times L^2(\mathbb{R}^d)$, $(W_T \otimes Z_S)(\phi_T, \phi_S)$ and $(W_T \otimes Z_S)(\psi_T, \psi_S)$ are centered Gaussian random variables with covariance

$$
\mathrm{Cov}[(W_T \otimes Z_S)(\phi_T, \phi_S), (W_T \otimes Z_S)(\psi_T, \psi_S)] = \int_0^T \phi_T \psi_T \, \mathrm{d}t \int_{\mathbb{R}^d} Z_S(\phi_S) Z_S(\psi_S) \, \mathrm{d}s.
$$

The space–time forcing term $W_T \otimes Z_S$ can be identified with a cylindrical Wiener process $\{\widetilde{W}_t\}_{t \in [0,T]}$ in $L^2(\mathbb{R}^d)$ through

$$
\widetilde{W}_t(\phi_S) = (W_T \otimes Z_S)(\mathbb{1}_{[0,t]}, \phi_S), \quad \phi_S \in L^2(\mathbb{R}^d), \quad t \in [0, T],
$$

where $\mathbb{1}_{[0,t]}$ is the indicator function over $[0,t]$. Moreover, we can write

$$
(W_T \otimes Z_S)(\phi_T, \phi_S) = \int_{\mathbb{R}^d} \left( \int_0^T \phi_T \, d\widetilde{W}_t \right) \phi_S \, ds \quad \phi_S \in L^2(\mathbb{R}^d), \quad t \in [0, T],
$$

which leads to the interpretation of $W_T \otimes Z_S$ as the time derivative of the cylindrical Wiener process $\{\widetilde{W}_t\}_{t \in [0,T]}$. This analogy allows
us to rewrite the SPDE (A.1) in the SDE form in infinite dimensions (Da Prato and Zabczyk, 1992)

$$
dX = -\frac{1}{c} \left[ (\kappa^2 - \nabla \cdot \mathbf{H} \nabla)^\alpha + \gamma \cdot \nabla \right] X dt + \frac{\tau}{\sqrt{c}} d\widetilde{W}_t.
$$

**Appendix B. Proof of Proposition 1**

We present here the proof of Proposition 1.

**Proof.** The covariance function of the spatial trace between $X(t, s)$ and $X(t, s')$ for a spatial lag $\mathbf{h} = s - s'$ does not depend on the
imaginary part of the spatial symbol function (Carrizo-Vergara et al., 2022), hence it can be written as

$$
\begin{align*}
\mathrm{Cov}(0, \mathbf{h}) &= \int_{\mathbb{R}^d} \int_{\mathbb{R}} \exp(i \mathbf{h} \cdot \xi) S(\xi, \omega) \, d\omega \, d\xi \\
&= \int_{\mathbb{R}^d} \exp(i \mathbf{h} \cdot \xi) \left[ \int_{\mathbb{R}} S(\xi, \omega) \, d\omega \right] \, d\xi \\
&= \int_{\mathbb{R}^d} \exp(i \mathbf{h} \cdot \xi) S_S(\xi) \, d\xi, \tag{B.3}
\end{align*}
$$

where $S(\xi, \omega)$ is the spectral density defined as

$$
S(\xi, \omega) = \frac{\tau^2}{(2\pi)^{d+1} [\omega^2 + c^{-2}(\kappa^2 + \xi^\top \mathbf{H} \xi)^{2\alpha}] c(\kappa^2 + \xi^\top \mathbf{H} \xi)^{\alpha_S}}.
$$

Integrating over $\omega$, we obtain the spatial spectral density

$$
\begin{align*}
S_S(\xi) &= \frac{\tau^2}{(2\pi)^d c (\kappa^2 + \xi^\top \mathbf{H} \xi)^{\alpha_S}} \int_{\mathbb{R}} \frac{1}{2\pi [\omega^2 + c^{-2}(\kappa^2 + \xi^\top \mathbf{H} \xi)^{2\alpha}]} \, d\omega \\
&= \frac{\tau^2}{(2\pi)^d c (\kappa^2 + \xi^\top \mathbf{H} \xi)^{\alpha_S}} \frac{1}{2 [c^{-2}(\kappa^2 + \xi^\top \mathbf{H} \xi)^{2\alpha}]^{1/2}} \\
&= \frac{\tau^2}{2(2\pi)^d (\kappa^2 + \xi^\top \mathbf{H} \xi)^{\alpha_{\text{tot}}}}. \tag{B.4}
\end{align*}
$$

Using the change of variable $\xi = \kappa \mathbf{H}^{-1/2} \mathbf{w}$ and plugging Eq. (B.4) into (B.3), we obtain

$$
\mathrm{Cov}(0, \mathbf{h}) = \frac{\tau^2}{2} \int_{\mathbb{R}^d} \frac{e^{i\mathbf{h}\cdot\xi}}{(2\pi)^d (\kappa^2 + \xi^\top \mathbf{H} \xi)^{\alpha_{\text{tot}}}} \mathrm{d}\xi
$$
$$
\begin{align*}
&= \frac{\tau^2}{2} \int_{\mathbb{R}^d} \frac{e^{i\mathbf{h}\cdot\mathbf{k} \mathbf{H}^{-1/2} \mathbf{w}} |\mathbf{k} \mathbf{H}^{-1/2}|}{(2\pi)^d (\kappa^2 + \kappa^2 \mathbf{w}^\top \mathbf{w})^{\alpha_{\text{tot}}}} \mathrm{d}\mathbf{w} \\
&= \frac{\tau^2}{2\kappa^2 (\alpha_{\text{tot}} - d/2) |\mathbf{H}|^{1/2}} \int_{\mathbb{R}^d} \frac{e^{i\mathbf{h}\cdot\mathbf{k} \mathbf{H}^{-1/2} \mathbf{w}}}{(2\pi)^d (1 + \mathbf{w}^\top \mathbf{w})^{\alpha_{\text{tot}}}} \mathrm{d}\mathbf{w} \\
&= \frac{\tau^2 \Gamma(\alpha_{\text{tot}} - d/2)}{2\Gamma(\alpha_{\text{tot}})(4\pi)^{d/2} \kappa^2 (\alpha_{\text{tot}} - d/2) |\mathbf{H}|^{1/2}} C_{\alpha_{\text{tot}} - d/2}^M (\kappa \|\mathbf{H}^{-1/2} \mathbf{h}\|).
\end{align*}
$$

The last result comes from the computation of $\int_{\mathbb{R}^d} (1 + \mathbf{w}^\top \mathbf{w})^{-\alpha_{\text{tot}}} \mathrm{d}\mathbf{w}$ with polar coordinates.

## Appendix C. Discretization of spatio-temporal advection–diffusion SPDE

We detail here the discretization scheme of the advection–diffusion spatio-temporal SPDE (5).

For the sake of a clearer exposition, we set $\mathbf{H} = \mathbf{I}$, $\alpha = 1$ and we consider a spatio-temporal white noise $Z(t, s) = W(t, s)$. The proof for the general case follows exactly the same lines as the proof below. The considered SPDE is

$$ \left[ \frac{\partial}{\partial t} + \frac{1}{c}(\kappa^2 - \Delta) + \frac{1}{c}\gamma \cdot \nabla \right] X(t, s) = \frac{\tau}{\sqrt{c}} W(t, s). \qquad (C.5) $$

For the discretization of the temporal derivative in Eq. (C.5), we opt for the implicit Euler scheme, which considers the differential equation

$$ \frac{\partial X(t)}{\partial t} = f(t, X), $$

with initial value $X^{(0)} = X(t_0)$. The method produces a sequence $\{X^{(k)}\}_{k=0}^{N_T}$, such that $X^{(k)}$ approximates $X(t_0 + kdt)$, where $dt$ is the time step size. The approximation reads

$$ X^{(k+1)} = X^{(k)} + d_t f(t^{(k+1)}, X^{(k+1)}). \qquad (C.6) $$

In the specific case of Eq. (C.5), the implicit Euler discretization step reads

$$ X^{(k+1)}(s) - X^{(k)}(s) + dt \left[ \frac{1}{c}(\kappa^2 - \Delta) + \frac{1}{c}\gamma \cdot \nabla \right] X^{(k+1)}(s) = \frac{\sqrt{dt\tau}}{\sqrt{c}} W_S^{(k+1)}(s), \qquad (C.7) $$

where $W_S^{(k+1)}(s)$ is a spatial white noise obtained by integrating out the temporal white noise.

For ease of notation, we denote $X^{(k+1)} = X^{(k+1)}(s)$, $X^{(k)} = X^{(k)}(s)$ and $W_S = W_S^{(k+1)}(s)$, since the spatial noise is independent of the temporal step $k$.

At each time step of the temporal discretization, a spatial finite element method is applied. In our case, we use the continuous Galerkin with Neumann boundary condition. The weak form of Eq. (C.7) is

$$ \begin{align}
\int_{\Omega} X^{(k+1)} v \, ds + \frac{dt}{c} \left( \int_{\Omega} \kappa^2 X^{(k+1)} v \, ds - \int_{\Omega} \Delta X^{(k+1)} v \, ds + \int_{\Omega} \gamma \cdot \nabla X^{(k+1)} v \, ds \right) &= \nonumber \\
&= \int_{\Omega} X^{(k)} v \, ds + \frac{\sqrt{dt\tau}}{\sqrt{c}} \int_{\Omega} v W_S(d\mathbf{s}), \quad \forall v \in \mathcal{V}, \tag{C.8}
\end{align} $$

where $\mathcal{V}$ is the Hilbert space in which we search the solution and $W_S(v)$ is the white noise applied to the test function $v$.

By applying Green's first identity, i.e., by writing

$$ \int_{\Omega} \Delta X^{(k+1)} v \, ds = - \int_{\Omega} \nabla X^{(k+1)} \cdot \nabla v \, ds + \int_{\partial \Omega} v \cdot (\nabla X^{(k+1)} \cdot \hat{\mathbf{n}}) \, d\sigma, $$

with $\hat{\mathbf{n}}$ being the normal vector on the boundary, and by simplifying the second term thanks to the Neumann boundary condition, we obtain

$$ \underbrace{\int_{\Omega} X^{(k+1)} v \, ds + \frac{dt}{c} \left( \int_{\Omega} \kappa^2 X^{(k+1)} v \, ds + \int_{\Omega} \nabla X^{(k+1)} \cdot \nabla v \, ds + \int_{\Omega} \gamma \cdot \nabla X^{(k+1)} v \, ds \right)}_{A(X^{(k+1)}, v)} \\
= \underbrace{\int_{\Omega} X^{(k)} v \, ds}_{C(X^{(k)}, v)} + \underbrace{\frac{\sqrt{dt\tau}}{\sqrt{c}} \int_{\Omega} v W_S(d\mathbf{s})}_{\mathcal{E}(v)}, \quad \forall v \in \mathcal{V}. $$

Let $\mathcal{V}_h$ be the space of finite element solutions spanned by the basis functions $\{\psi_i\}_{i=1}^{N_S}$. The Galerkin method allows us to find an approximated solution $X_h^{(k+1)} \in \mathcal{V}_h \subset \mathcal{V}$ to the SPDE, such that

$$ \mathcal{A}(X_h^{(k+1)}, v_h) = C(X_h^{(k)}, v_h) + \mathcal{E}(v_h) \quad \forall v_h \in \mathcal{V}_h. \qquad (C.9) $$
The functions $X_h^{(k+1)}$, $X_h^{(k)}$ and $v_h$ are linear combinations of the basis functions, with

$$X_h^{(k+1)} = \sum_{i=1}^{N_S} x_i^{(k+1)} \psi_i; \quad X_h^{(k)} = \sum_{i=1}^{N_S} x_i^{(k)} \psi_i; \quad v_h = \sum_{i=1}^{N_S} v_i \psi_i.$$

Because of the linearity in the first argument of $\mathcal{A}(\cdot, \cdot)$ and $\mathcal{C}(\cdot, \cdot)$, we get

$$\sum_{i=1}^{N_S} \mathcal{A}(\psi_i, v_h) x_i^{(k+1)} = \sum_{i=1}^{N_S} \mathcal{C}(\psi_i, v_h) x_i^{(k)} + \mathcal{E}(v_h), \quad \forall v_h \in \mathcal{V}_h, \qquad (C.10)$$

where

$$\mathcal{A}(\psi_i, v_h) = \mathcal{M}(\psi_i, v_h) + \frac{dt}{c} (\mathcal{K}(\psi_i, v_h) + \mathcal{B}(\psi_i, v_h))$$

$$\mathcal{C}(\psi_i, v_h) = \mathcal{M}(\psi_i, v_h),$$

with $\mathcal{K}(\psi_i, v_h) = \kappa^2 \mathcal{M}(\psi_i, v_h) + \mathcal{G}(\psi_i, v_h)$. Here, $\mathcal{M}$ and $\mathcal{G}$ are the mass and stiffness operators, respectively $\mathcal{M}(v, w) = \int_{\Omega} vw ds$ and $\mathcal{G}(v, w) = \int_{\Omega} \nabla v \cdot \nabla w ds$. $B$ is the advection operator, i.e., $\mathcal{B}(v, w) = \int_{\Omega} \gamma \cdot \nabla vw ds$. Finally, $\mathcal{E}$ is the operator of the form $\mathcal{E}(v) = \frac{\sqrt{dt}}{\sqrt{c}} \int_{\Omega} v W_S ds$.

Since any $v_h$ can be written as a linear combination of basis functions, the formulation (C.10) is equivalent to

$$\sum_{i=1}^{N_S} \mathcal{A}(\psi_i, \psi_j) x_i^{(k+1)} = \sum_{i=1}^{N_S} \mathcal{C}(\psi_i, \psi_j) x_i^{(k)} + \mathcal{E}(\psi_j), \quad \forall j. \qquad (C.11)$$

We define $\mathbf{M} = [M_{ij}]_{i,j=1}^{N_S} = [\mathcal{M}(\psi_i, \psi_j)]_{i,j=1}^{N_S}$, $\mathbf{G} = [G_{ij}]_{i,j=1}^{N_S} = [\mathcal{G}(\psi_i, \psi_j)]_{i,j=1}^{N_S}$, $\mathbf{B} = [B_{ij}]_{i,j=1}^{N_S} = [\mathcal{B}(\psi_i, \psi_j)]_{i,j=1}^{N_S}$ the mass, stiffness and advection matrices, respectively.

$\mathcal{E}(\psi_j)$ is a Gaussian random variable with expectation 0 and covariance equal to

$$\begin{aligned} \mathrm{Cov}(\mathcal{E}(\psi_i), \mathcal{E}(\psi_j)) &= \frac{dt}{c} \mathrm{Cov}[W_S(\psi_i), W_S(\psi_j)] \\ &= \frac{dt}{c} \int_{\Omega} \psi_i \psi_j \, ds = \frac{dt}{c} M_{ij}, \end{aligned}$$

by following the definition of white noise in Eq. (2).

If $z^{(k+1)}$ is a $(N_S)$-Gaussian vector such that $z^{(k+1)} \sim \mathcal{N}(0, \mathbf{I}_{N_S})$, $\mathbf{x}^{(k+1)}$ is the vector containing the values $\{x_i^{(k+1)}\}_{i=1}^{N_S}$ and $\mathbf{x}^{(k)}$ is the vector containing the values $\{x_i^{(k)}\}_{i=1}^{N_S}$, then the sparse linear system corresponding to Eq. (C.11) reads

$$\mathbf{M} \mathbf{x}^{(k+1)} + \frac{dt}{c} (\mathbf{K} + \mathbf{B}) \mathbf{x}^{(k+1)} = \mathbf{M} \mathbf{x}^{(k)} + \frac{\sqrt{dt}}{\sqrt{c}} \mathbf{M}^{1/2} \mathbf{z}^{(k+1)}, \qquad (C.12)$$

where $\mathbf{K} = \kappa^2 \mathbf{M} + \mathbf{G}$ and $\mathbf{M}^{1/2}$ is any matrix such that $\mathbf{M}^{1/2} \mathbf{M}^{1/2} = \mathbf{M}$.

When the spatial noise is colored, i.e. $Z_S(s)$, the right-hand side operator $\mathcal{E}(v)$ becomes $\mathcal{E}_S(v)$, defined as

$$\mathcal{E}_S(v) = \frac{\sqrt{dt}}{\sqrt{c}} \int_{\Omega} v Z_S ds,$$

and it satisfies

$$\mathcal{E}_S(v_h) = \sum_{i=1}^{N_S} \mathcal{M}(\psi_i, v_h) z_{S,i}.$$

Hence,

$$\sum_{i=1}^{N_S} \mathcal{A}(\psi_i, \psi_j) x_i^{(k+1)} = \sum_{i=1}^{N_S} \mathcal{C}(\psi_i, \psi_j) x_i^{(k)} + \sum_{i=1}^{N_S} \mathcal{M}(\psi_i, \psi_j) z_{S,i}, \quad \forall j.$$

If $z_S = \{z_{S,i}\}_{i=1}^{N_S}$ has precision matrix equal to $Q_S$, then the sparse linear system is

$$\mathbf{M} \mathbf{x}^{(k+1)} + \frac{dt}{c} (\mathbf{K} + \mathbf{B}) \mathbf{x}^{(k+1)} = \mathbf{M} \mathbf{x}^{(k)} + \frac{\sqrt{dt}}{\sqrt{c}} \mathbf{M} \mathbf{L}_S^\top \mathbf{z}^{(k+1)}, \qquad (C.13)$$

where $\mathbf{z}^{(k+1)} \sim \mathcal{N}(0, \mathbf{I}_{N_S})$ and $\mathbf{L}_S$ is the Cholesky decomposition of $\mathbf{Q}_S^{-1}$.

**Remark 3.** When the diffusion term includes an anisotropy matrix $\mathbf{H}$, i.e., when $\Delta$ is replaced by $\nabla \cdot \mathbf{H} \nabla$, the stiffness operator becomes $\mathcal{G}(v, w) = \int_{\Omega} \mathbf{H} \nabla v \cdot \nabla w \, ds$, and the stiffness matrix changes consequently.
## Appendix D. Proof of Proposition 6

We present here the proof of Proposition 6.

**Proof.** Let us denote $\mathbf{x}_{0:N_T} = [\mathbf{x}^{(0)}, \dots, \mathbf{x}^{(N_T)}]^\mathrm{T}$ the vector containing all spatial solutions until time step $N_T$. Then,

$$ \mathbf{x}_{0:N_T} = \mathbf{R} \left( \mathbf{z}_{1:N_T}^{(0)} \right), $$

with $\mathbf{z}_{1:N_T} = [\mathbf{z}^{(1)}, \dots, \mathbf{z}^{(N_T)}]^\mathrm{T}$ and

$$ \mathbf{R} = \begin{pmatrix} \mathbf{I}_{N_S} & 0 & 0 & \dots & 0 \\ \mathbf{D} & \mathbf{E} & 0 & 0 & 0 \\ \mathbf{D}^2 & \mathbf{D} \mathbf{E} & \mathbf{E} & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \vdots & \vdots & \vdots & \ddots & 0 \\ \vdots & \vdots & \mathbf{D}^2 & \mathbf{D} & \mathbf{E} \end{pmatrix}. $$

$\mathbf{R}$ has a block structure which allows easy computation of its inverse

$$ \mathbf{R}^{-1} = \begin{pmatrix} \mathbf{I}_{N_S} & 0 & 0 & \dots & 0 \\ -\mathbf{E}^{-1} \mathbf{D} & \mathbf{E}^{-1} & 0 & 0 & 0 \\ 0 & -\mathbf{E}^{-1} \mathbf{D} & \mathbf{E}^{-1} & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & \dots & -\mathbf{E}^{-1} \mathbf{D} & \mathbf{E}^{-1} \end{pmatrix}. $$

The precision matrix of $\mathbf{x}_{0:N_T}$ is thus

$$ \mathbf{Q} = \mathbf{R}^{-1} \begin{pmatrix} \Sigma^{-1} & 0 & \dots & 0 \\ 0 & \mathbf{I}_{N_S} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \mathbf{I}_{N_S} \end{pmatrix} \mathbf{R}^{-1}. $$

By denoting $\mathbf{F} = \mathbf{E}\mathbf{E}^\mathrm{T}$, the global precision matrix reads

$$ \mathbf{Q} = \begin{pmatrix} \Sigma^{-1} + \mathbf{D}^\mathrm{T} \mathbf{F}^{-1} \mathbf{D} & -\mathbf{D}^\mathrm{T} \mathbf{F}^{-1} & 0 & \dots & 0 \\ -\mathbf{F}^{-1} \mathbf{D} & \mathbf{F}^{-1} + \mathbf{D}^\mathrm{T} \mathbf{F}^{-1} \mathbf{D} & -\mathbf{D}^\mathrm{T} \mathbf{F}^{-1} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & -\mathbf{F}^{-1} \mathbf{D} & \mathbf{F}^{-1} + \mathbf{D}^\mathrm{T} \mathbf{F}^{-1} \mathbf{D} & -\mathbf{D}^\mathrm{T} \mathbf{F}^{-1} & 0 \\ 0 & 0 & -\mathbf{F}^{-1} \mathbf{D} & \mathbf{F}^{-1} & 0 \end{pmatrix}. $$

By replacing the values of $\mathbf{D}$ and $\mathbf{F}$ and by defining $\mathbf{J} = [\mathbf{M} + \frac{dt}{c}(\mathbf{K} + \mathbf{B} + \mathbf{S})]$, we obtain

$$ \mathbf{Q} = \frac{c}{r^2 dt} \begin{pmatrix} \Sigma^{-1} + Q_S & -Q_S M^{-1} J & 0 & \dots & 0 \\ -J^\mathrm{T} M^{-1} Q_S & J^\mathrm{T} M^{-1} Q_S M^{-1} J + Q_S & -Q_S M^{-1} J & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & -J^\mathrm{T} M^{-1} Q_S & J^\mathrm{T} M^{-1} Q_S M^{-1} J + Q_S & -Q_S M^{-1} J & 0 \\ 0 & 0 & -J^\mathrm{T} M^{-1} Q_S & J^\mathrm{T} M^{-1} Q_S M^{-1} J & 0 \end{pmatrix}. $$

## References

Ailliot, P., Baxevani, A., Cuzol, A., Monbet, V., Raillard, N., 2011. Space-time models for moving fields with an application to significant wave height fields. Environmetrics 22 (3), 354–369.

Allard, D., Clarotto, L., Emery, X., 2022. Fully nonseparable Gneiting covariance functions for multivariate space-time data. Spatial Stat. 52, 100706.

Bakka, H., 2022. How to solve the stochastic partial differential equation that gives a matérn random field using the finite element method. working paper or preprint: arXiv:1803.03765.

Banerjee, S., Carlin, B.P., Gelfand, A.E., 2014. Hierarchical Modeling and Analysis for Spatial Data, second ed. Chapman & Hall/CRC Press. Boca Raton, FL.

Bank, R.E., Bürgler, J.F., Fichtner, W., Smith, R.K., 1990. Some upwinding techniques for finite element approximations of convection-diffusion equations. Numer. Math. 58 (1), 185–202.

Benoit, L., Allard, D., Mariethoz, G., 2018. Stochastic rainfall modeling at sub-kilometer scale. Water Resour. Res. 54 (6), 4108–4130.

Beyer, H.G., Costanzo, C., Heinemann, D., 1996. Modifications of the Heliosat procedure for irradiance estimates from satellite images. Sol. Energy 56 (3), 207–212.

Bolin, D., Kirchner, K., 2020. The Rational SPDE Approach for Gaussian Random Fields With General Smoothness. J. Comput. Graph. Statist. 29 (2), 274–285.

Bolin, D., Simas, A.B., Xiong, Z., 2024. Covariance-based rational approximations of fractional SPDEs for computationally efficient Bayesian inference. J. Comput. Graph. Statist. 33 (1), 64–74.

Bourrette, M., Allard, D., Porcu, E., 2016. A flexible class of non-separable cross-covariance functions for multivariate space-time data. Spatial Stat. 18, 125–146.

Cameletti, M., Lindgren, F., Simpson, D., Rue, H., 2011. Spatio-temporal modeling of particulate matter concentration through the SPDE approach. AStA Adv. Stat. Anal. 97.

Carrizo-Vergara, R., Allard, D., Desassis, N., 2022. A general framework for SPDE-based stationary random fields. Bernoulli 28 (1), 1–32.

Chebyshev, P.L., 1853. Théorie des mécanismes connus sous le nom de parallélogrammes. Imprimerie de l'Académie impériale des sciences, St. Pétersbourg.

Chen, W., Genton, M.G., Sun, Y., 2021. Space-Time Covariance Structures and Models. Annu. Rev. Stat. Appl. 8 (1), 191–215.

Chilès, J.-P., Delfiner, P., 1999. Geostatistics: Modeling Spatial Uncertainty. John Wiley & Sons, New York.
Cressie, N., Huang, H.-C., 1999. Classes of nonseparable, spatio-temporal stationary covariance functions. J. Amer. Statist. Assoc. 94 (448), 1330-1340.

Cressie, N., Wikle, C.K., 2011. Statistics for Spatio-Temporal Data. Wiley.

Da Prato, G., Zabczyk, J., 1992. Stochastic Equations in Infinite Dimensions. Encyclopedia of Mathematics and its Applications, Cambridge University Press.

Erisman, A.M., Tinney, W.F., 1975. On computing certain elements of the inverse of a sparse matrix. Commun. ACM 18 (3), 177-179.

Fuglstad, G.-A., Simpson, D., Lindgren, F., Rue, H., 2015. Does non-stationary spatial data always require non-stationary random fields? Spatial Stat. 14 (Part C), 505-531.

Genton, M., 2001. Classes of kernels for machine learning: A statistics perspective. J. Mach. Learn. Res. 2, 299-312.

Gneiting, T., 2002. Nonseparable, stationary covariance functions for space-time data. J. Amer. Statist. Assoc. 97 (458), 590-600.

Gschwind, B., Wald, L., Blanc, P., Lefèvre, M., Schroedter-Homscheidt, M., Arola, A., 2019. Improving the McClear model estimating the downwelling solar radiation at ground level in cloud-free conditions - McClear-v3. Meteorol. Z. 28 (2), 147-163.

Heaton, M., Datta, A., Finley, A., Furrer, R., Guinness, J., Guhaniyogi, R., Gerber, F., Gramacy, R., Hammerling, D., Katzfuss, M., Lindgren, F., 2019. A case study competition among methods for analyzing large spatial data. J. Agric. Biol. Environ. Stat. 24, 398-425.

Huang, H., Castruccio, S., Genton, M.G., 2022. Forecasting high-frequency spatio-temporal wind power with dimensionally reduced echo state networks. J. R. Stat. Soc. Ser. C. Appl. Stat. 71 (2), 449-466.

Hughes, T., Brooks, A., 1979. A multidimensional upwind scheme with no cross-wind diffusion in finite element methods for advection dominated flows. p. 34, ASME Publication ASME-AMD.

Hutchinson, M., 1990. A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Comm. Statist. Simulation Comput. 19 (2), 433-450.

Krainski, E., Gómez Rubio, V., Bakka, H., Lenzi, A., Castro-Camilo, D., Simpson, D., Lindgren, F., Rue, H., 2018. Advanced Spatial Modeling with Stochastic Partial Differential Equations Using R and INLA. Chapman & Hall/CRC Press, Boca Raton, FL.

Kreuwel, F.P., Knap, W.H., Visser, L.R., van Sark, W.G., de Arellano, J.V.-G., van Heerwaarden, C.C., 2020. Analysis of high frequency photovoltaic solar energy fluctuations. Sol. Energy 206, 381-389.

Lenzi, A., Genton, M.G., 2020. Spatiotemporal probabilistic wind vector forecasting over Saudi Arabia. Ann. Appl. Stat. 14 (3), 1359-1378.

Lindgren, F., Bakka, H., Bolin, D., Krainski, E., Rue, H., 2023. A diffusion-based spatio-temporal extension of Gaussian Matérn fields. working paper or preprint: arXiv:2006.04917.

Lindgren, F., Bolin, D., Rue, H., 2022. The SPDE approach for Gaussian and non-Gaussian fields: 10 years and still running. Spatial Stat. 50, 100599, Special Issue: The Impact of Spatial Statistics.

Lindgren, F., Rue, H., 2015. Bayesian Spatial Modelling with R-INLA. J. Stat. Softw. 63 (19), 1-25.

Lindgren, F., Rue, H.a., Lindström, J., 2011. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. J. R. Stat. Soc. B Stat. Methodol. 73 (4), 423-498.

Liu, X., Yeo, K., Lu, S., 2020. Statistical modeling for spatio-temporal data from stochastic convection-diffusion processes. J. Amer. Statist. Assoc. 1-37.

Macke, A., Seifert, P., Baars, H., Barthlott, C., Beekmans, C., Behrendt, A., Bohn, B., Brueck, M., Bühl, J., Crewell, S., Damian, T., Deneke, H., Düsing, S., Foth, A., Di Girolamo, P., Hammann, E., Heinze, R., Hirsikko, A., Kalisch, J., Kalthoff, N., Kinne, S., Kohler, M., Löhnert, U., Madhavan, B.L., Maurer, V., Muppa, S.K., Schween, J., Serikov, I., Siebert, H., Simmer, C., Späth, F., Steinke, S., Träummer, K., Trömel, S., Wehner, B., Wieser, A., Wulfmeyer, V., Xie, X., 2017. The hd(CP)² observational prototype experiment (HOPE) – an overview. Atmos. Chem. Phys. 17 (7), 4887-4914.

Martínez-Hernández, I., Genton, M.G., 2023. Surface time series models for large spatio-temporal datasets. Spatial Statistics 53, 100718.

Mekuria, G., Rao, J., 2016. Adaptive finite element method for steady convection-diffusion equation. Am. J. Comput. Math. 06, 275-285.

Nocedal, J., Wright, S.J., 2006. Numerical Optimization, 2e ed. Springer, New York, NY, USA.

Okasaki, C., Hooten, M.B., Berdahl, A.M., 2022. Source reconstruction for spatio-temporal physical statistical models. Spatial Stat. 52, 100707.

Oumbe, A., Qu, Z., Blanc, P., Lefèvre, M., Wald, L., Cros, S., 2014. Decoupling the effects of clear atmosphere and clouds to simplify calculations of the broadband solar irradiance at ground level. Geosci. Model Dev. 7 (4), 1661-1669.

Paciorek, C.J., Yanosky, J.D., Puett, R.C., Laden, F., Suh, H.H., 2009. Practical large-scale spatio-temporal modeling of particulate matter concentrations. Ann. Appl. Stat. 3 (1), 370-397.

Pereira, M., 2023. A note on spatio-temporal random fields on meshed surfaces defined from advection-diffusion SPDEs. working paper or preprint: URL https://hal.science/hal-04132148.

Pereira, M., Desassis, N., 2019. Efficient simulation of Gaussian Markov random fields by Chebyshev polynomial approximation. Spatial Stat. 31, 100359.

Pereira, M., Desassis, N., Allard, D., 2022. Geostatistics for large datasets on Riemannian manifolds: A matrix-free approach. J. Data Sci. 20 (4), 512-532.

Porcu, E., Furrer, R., Nychka, D., 2021. 30 Years of space-time covariance functions. Wiley Interdiscip. Rev. Comput. Stat. 13 (2), e1512.

Porcu, E., Gregori, P., Mateu, J., 2006. Nonseparable stationary anisotropic space-time covariance functions. Stoch. Environ. Res. Risk Assess. 21, 113-122.

Powell, P.D., 2011. Calculating Determinants of Block Matrices. working paper or preprint: arXiv:1112.4379.

Quarteroni, A., 2008. Modellistica Numerica per Problemi Differenziali. Springer Milano.

Roques, L., Allard, D., Soubeyrand, S., 2022. Spatial statistics and stochastic partial differential equations: A mechanistic viewpoint. Spatial Stat. 50, 100591.

Rue, H., Held, L., 2005. Gaussian Markov random fields: theory and applications, first ed. In: Monographs on statistics and applied probability, vol. 104, Chapman & Hall/CRC.

Rue, H., Martino, S., Chopin, N., 2009. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. J. R. Stat. Soc. Ser. B Stat. Methodol. 71 (2), 319-392.

Salvaña, M.L.O., Genton, M.G., 2021. Lagrangian Spatio-Temporal Nonstationary Covariance Functions. In: Advances in Contemporary Statistics and Econometrics: Festschrift in Honor of Christine Thomas-Agnan. Springer International Publishing, Cham, pp. 427-447.

Särkkä, S., Solin, A., Hartikainen, J., 2013. Spatiotemporal Learning via Infinite-Dimensional Bayesian Filtering and Smoothing: A Look at Gaussian Process Regression Through Kalman Filtering. IEEE Signal Process. Mag. 30, 51-61.

Schade, N.H., Macke, A., Sandmann, H., Stick, C., 2007. Enhanced solar global irradiance during cloudy sky conditions. Meteorol. Z. 16, 295-303.

Sharrock, L., Kantas, N., 2022. Joint online parameter estimation and optimal sensor placement for the partially observed stochastic advection-diffusion equation. SIAM/ASA J. Uncertainty Quant. 10, 55-95.

Sigrist, F., Künsch, H.R., Stahel, W.A., 2011. An autoregressive spatio-temporal precipitation model. Procedia Environ. Sci. 3, 2-7, 1st Conference on Spatial Statistics 2011 – Mapping Global Change.

Sigrist, F., Künsch, H.R., Stahel, W.A., 2012. A dynamic nonstationary spatio-temporal model for short term prediction of precipitation. Ann. Appl. Stat. 6 (4), 1452-1477.

Sigrist, F., Künsch, H.R., Stahel, W.A., 2015. Stochastic partial differential equation based modelling of large space-time data sets. J. R. Stat. Soc. Ser. B Stat. Methodol. 77 (1), 3-33.

Stein, M.L., 2005. Space-Time Covariance Functions. J. Amer. Statist. Assoc. 100 (469), 310-321.

Takahashi, K., Fagan, J., Chin, M., 1973. Formation of a sparse bus impedance matrix and its application to short circuit study. In: 8th PICA Conf. Proc. June 4-6, Minneapolis, Minn.

Vabishchevich, P.N., 2015. Numerically solving an equation for fractional powers of elliptic operators. J. Comput. Phys. 282, 289-302.

Whittle, P., 1954. On stationary processes in the plane. Biometrika 41 (3-4), 434-449.

Whittle, P., 1963. Stochastic processes in several dimensions. Bull. Inst. Int. Stat. 40, 974-994.

Wikle, C.K., Cressie, N., 1999. A dimension-reduced approach to space-time Kalman filtering. Biometrika 86 (4), 815-829.

Wikle, C., Hooten, M., 2010. A general science-based framework for dynamical spatio-temporal models. TEST 19, 417-451.

Young, D.M., 1971. Iterative Solution of Large Linear Systems. Academic Press.