# Boundary Conditions for the Parametric Kalman Filter forecast
version of: ,

M. Sabathier
ONERA, Toulouse, France
<EMAIL_ADDRESS>
&O. Pannekoucke
INPT-ENM, Toulouse, France
CNRM, Université de Toulouse, Météo-France, CNRS, Toulouse, France
CERFACS, Toulouse, France
<EMAIL_ADDRESS>
&V. Maget
ONERA, Toulouse, France
<EMAIL_ADDRESS>
&N. Dahmen
ONERA, Toulouse, France
<EMAIL_ADDRESS>
Corresponding author

###### Abstract

This paper is a contribution to the exploration of the parametric Kalman
filter (PKF), which is an approximation of the Kalman filter, where the error
covariances are approximated by a covariance model. Here we focus on the
covariance model parameterized from the variance and the anisotropy of the
local correlations, and whose parameters dynamics provides a proxy for the
full error-covariance dynamics. For this covariance model, we aim to provide
the boundary condition to specify in the prediction of PKF for bounded
domains, focusing on Dirichlet and Neumann conditions when they are prescribed
for the physical dynamics. An ensemble validation is proposed for the
transport equation and for the heterogeneous diffusion equation over a bounded
1D domain. This ensemble validation requires to specify the auto-correlation
time-scale needed to populate boundary perturbation that leads to prescribed
uncertainty characteristics. The numerical simulations show that the PKF is
able to reproduce the uncertainty diagnosed from the ensemble of forecast
appropriately perturbed on the boundaries, which show the ability of the PKF
to handle boundaries in the prediction of the uncertainties. It results that
Dirichlet condition on the physical dynamics implies Dirichlet condition on
the variance and on the anisotropy.

_K_ eywords Data assimilation $\cdot$ Parametric Kalman filter $\cdot$
Boundary conditions

## Plain Language Summary

This work addresses the question of the uncertainty prediction in bounded
domains. It contributes to explore a theoretical formulation of the
uncertainty prediction, the parametric Kalman filter, that opens the way to
data assimilation in real applications where the boundaries are important,
such as in radiation belts predictions, air quality, or wild-land fire ; or
the exploration of the coupling of the uncertainty in the atmosphere-ocean
coupled system.

##  1 Introduction

Uncertainty prediction is a challenging topic, important in probabilistic
forecasting as well as in data assimilation (DA) e.g. in variational data
assimilation to improve the specification of the background term [21, 6]. One
of the main theoretical backbone in DA is given by the Kalman filter equations
that applies for linear dynamics, but that cannot be used for large systems
because the numerical cost to predict the error covariance becomes
prohibitive. Hence, approximations of the KF have been proposed, as the
ensemble Kalman filter (EnKF) where the error covariance matrices are
approximated from ensemble estimation [14].

The EnKF is widely used and has shown to perform well for many applications in
geosciences e.g. for the weather prediction [18], or for the radiation belts
prediction [4]. Radiation belts dynamics modeling consist in estimating
quantitatively the fluxes of high energetic electrons and protons trapped in
the Earth magnetic field using a typical advection-diffusion equation. This
region spans from 1 Earth Radius up to 8, thus encompassing all typical
satellites orbits, with which such particles can strongly interact and induce
minor to critical onboard anomalies. Compared with global prediction,
radiation belts predictions are performed on a limited and non-periodic domain
where the boundary imposes conditions to the dynamics of the electrons and
protons. Indeed, on one side, the outer boundary condition is considered as
the prime access for fresh materials, coming from the so-called magneto-tail,
and is typically modeled as an imposed Dirichlet condition at this altitude (8
Earth radii) that can evolve as a function of solar activity (e.g. energy
spectrum reshaping from time to time) [22]. On the other side, close to the
Earth, the atmosphere implies a necessary fixed Dirichlet condition too, as
all radiation belts particles coming down there are absorbed (e.g.
distribution always equal to 0). Finally, for low energy boundary we expect to
rely on a Neumann condition to limit naturally any escape of particles or
artificial source. Nonetheless, when performed on limited area models,
atmospheric prediction also present such kind of boundaries.

In weather forecasting over local-area domains, the importance of the boundary
has been identified from sensitivity studies showing that errors may be traced
back to the boundaries [13]. Hence, inaccurate boundary conditions may result
in large forecast errors [49, 50, 51]. However, in preliminary studies on the
potential of EnKF applied to local-area models (LAM), no perturbations of the
boundary conditions were considered [43, 12], leading to a loss of variance as
lead time increases [27, 28]. Note that, this loss of variance also occurs for
the Kalman filter prediction step e.g. in a transport equation over 1D
limited-area domain, when the boundary condition is not well formalized in the
dynamics, the inflow of variance is null at the inlet of the domain (see A).
This issue can be solved from the formalism of the Kalman filter as used in
optimal control, by considering that the boundary is an uncertain command as
done in e.g. heat conduction estimation [42]. In the community of the EnKF for
LAMs, the loss of variance is avoided by perturbing the boundary conditions
[28], e.g. by modeling the spatial and temporal covariance relationship of the
boundary condition based on the use of multivariate covariance model or by
using boundary conditions which derives from an ensemble on a larger domain
[47], even if for the latter approach, the consistency across multiple domains
is difficult to handle [19, sec. 6.a]. Note that in variational data
assimilation for LAMs, it is possible to include the uncertainty on the
boundary 4DVar through the adjunction of a cost function term that constraints
the increment or its tendency at the boundary [16], which is a way to avoid
the loss of variance above mentioned for the EnKF and for the KF.

In radiation belts where there is no larger model and thus no ensemble from
which perturbation on the boundary can be introduced, we are interested by
sampling boundary perturbations from a spatio-temporal covariance model.
However, because the solar activity is non-stationary, we need to develop
perturbations of the boundary conditions that are heterogeneous and non-
stationary i.e. with varying time auto-correlation, in place to the auto-
regressive model often used but for which the time auto-correlation is
constant [47]. Moreover, we want to reduce the numerical cost of using an EnKF
while keeping the uncertainty dynamics it provides. This motivation can be
also be encountered in the development of the assimilation scheme for
operational model in air quality, that deals with hundreds of chemical
species, and for which the numerical cost of an integration is large. For
instance, in the operational air quality model MOCAGE developed at Météo-
France, the assimilation relies on a 3D-FGAT approach [1, 24], without
ensemble, and it would be interesting to introduce flow-dependency of the
background term while accounting for the transport of the boundary uncertainty
during the assimilation at regional scale.

Recently another approximation for the KF, different from the EnKF, has been
introduced, the parametric Kalman filter (PKF), where the error covariance
matrices are approximated from a covariance model [38]. In the PKF, the
dynamics of the parameters provides a proxy for the dynamics of the full
covariance matrix. For instance, covariance model parameterized from the
variance and the anisotropy of the local correlation functions are able to
predict the dynamics of the covariance matrix for transport equations [9], but
at a numerical cost equivalent to three time the integration of the transport
[32, 35]. In addition, the PKF provides an understanding of the uncertainties
dynamics, which is less obvious to determine from an ensemble estimation
alone. For instance, when considering a transport, the PKF equations for the
error statistics consist in: the advection of the error variance without
source term, meaning that the magnitude of the error variance is conserved ;
and the advection of the anisotropy plus a source term that corresponds to the
deformation of the anisotropy due to the local shear [35]. As another example,
when considering now a conservative dynamics or a non-linear Burgers dynamics
(in the tangent-linear approximation for the error statistics), a source term
appears in the PKF equation for the error variance which is related to high
gradient area that can be interpreted as a production term similar to the one
encountered in turbulence [32]. These examples of physical interpretations of
the dynamics of the error statistics, as it is provided by the PKF equations,
are less easy to deduce from the ensemble estimation while the effects of
these processes can be observed. Note that data-driven modeling leveraged on
deep-learning, especially physics-informed neural network [40], can be
considered to determine the dynamics of the error statistics [37]. On the
other hand, it is not necessarily easy to interpret all terms in the PKF
equations e.g. when considering a diffusion equation, for which the PKF shows
a coupling between the error variance and the anisotropy, it is hard to
physically explain terms in square of the gradient of the anisotropy. In this
way, the PKF can contribute to improve the understanding of other topics of
interest in data assimilation e.g. in the characterization of the dynamics of
the model-error covariance due to the spatial and temporal discretization of a
transport equation [34, 25]. The PKF is interesting to provide a flow
dependency of the background error covariance matrix in variational data
assimilation e.g. by providing a proxy for the diffusion tensor of the day in
the covariance model based on the diffusion equation [38] ; or by providing
the variance of the day in presence of model error thanks to the PKF dynamics
of the variance, as in the assimilation of GOSAT methane observations using
the hemispherical CMAQ [48]. While the PKF has been mainly developed in
univariate statistics, an exploration has begun to develop a multivariate
statistics formulation [39]. At a more theoretical level, the parametric
formulation has been considered to describe the covariance dynamics of white
noise error in a conservative equation [15]. Hence, the PKF approach seems
promising to tackle the issues encountered in the development of the data
assimilation for radiation belts and air quality.

Until now, the PKF has been explored mainly on periodic 1D or 2D domains,
where it has been shown to reproduce interesting features of the uncertainty
dynamics in linear problem e.g. for the transport [35], as well as for non-
linear dynamics at the second order e.g. for the non-linear advection-
diffusion equation [32]. Note that the non-linear extension of the PKF relies
on a tangent-linear evolution of error, with a feedback of the uncertainty on
the dynamics of the mean through the error-variance (i.e. a fluctuation-mean
interaction), leading to a Gaussian second-order filter [20, Sec. 9.3].
However, to go ahead towards real applications, and especially applications in
bounded domains, appropriate specification of boundary conditions of the error
statistics is needed for the PKF dynamics. To do so, we propose to explore the
specification of the boundary conditions for the PKF when Dirichlet and
Neumann conditions are considered in the physical dynamics. This exploration
focuses on two dynamics of interest for our applications: the transport
equation e.g. for air quality or weather prediction ; and the diffusion
equation e.g. for radiation belts prediction or uncertainty dynamics in
boundary layer for air quality. The paper focuses on the forecast step, and
the assimilation step is not addressed here.

The paper is organized as follows. First, the background of the parametric
Kalman filter (PKF) is reminded in Section 2. Then, Section 3 details how to
specify the PKF conditions at the boundary for the forecast for the Dirichlet
and the Neumann conditions. The ensemble validation of the boundary conditions
for the PKF needs an ensemble of forecasts. To do so, an intermediate Section
4 will details how to specify boundary conditions in an EnKF experiment that
produces desired error statistics. This is an important contribution of the
paper in order to validate the specification of the boundary conditions of the
PKF, where the numerical validation is presented in Section 5. Conclusions and
perspectives are given in the last Section 6.

##  2 Background on the PKF forecast step

This section gives a self contained introduction to the parametric Kalman
filter, applied for a particular covariance model. First, the prediction step
of the Kalman filter applied on a linear dynamics is reminded. Then, the
formalism of the PKF is presented, followed by the illustration on two
dynamics: the transport equation, important in geosciences, and the diffusion
equation important in radiation belt dynamics community.

###  2.1 Kalman filter forecast step

Here we consider the prediction of a univariate physical field
$\chi(t,\mathbf{x})$ defined on a domain $\Omega$ of dimension $d$ and
coordinate system $\mathbf{x}=(x^{i})_{i\in[1,d]}$, whose dynamics is given by

$\partial_{t}\chi=\mathcal{M}(\partial\chi),$ (1)

where $\mathcal{M}$ stands for a function of the state $\chi$ and of its
spatial derivatives, $\partial\chi$, which is a shorthand for the partial
derivative with respect to the spatial coordinates at any arbitrary orders,
with the convention that order zero denotes the field $\chi$ itself.
Thereafter, for the sake of simplicity, $\mathcal{M}$ is assumed linear but
the formalism extends to the non-linear framework [32, 36]. Note that $\chi$
can be either continuous or discrete (the discretized version of the
continuous field): the discrete case leads to matrix algebra relations e.g.
$\mathcal{M}$ is replaced by its matrix formulation $\mathbf{M}$.

Hence, for a time integration over a time window $[0,T]$, the dynamics Eq. (1)
is written as

$\chi(T)=\mathbf{M}_{T\leftarrow 0}\chi(0),$ (2)

where $\mathbf{M}_{T\leftarrow 0}$ denotes the propagator associated with the
time integration of Eq. (1) over $[0,T]$. Compared to the classical
formulation of the Kalman filter equations, we do not consider here the model
error because it will be made negligible in the numerical experiments
performed in this study by considering a sufficiently fine space-time
resolution.

In real applications, the spatio-temporal heterogeneity of the observation
network, as well as the model error, imply that $\chi$ is not known exactly.
The true state of the system is denoted by $\chi^{t}$. The analysis state,
that is the estimation of the true state knowing the observations until a
given time, is denoted by $\chi^{a}$. The deviation of the analysis state from
the truth is the analysis error, $e^{a}=\chi^{a}-\chi^{t}$, and is often
modeled as a random Gaussian vector of zero mean and covariance matrix
$\mathbf{P}^{a}=\mathbb{E}\left[e^{a}(e^{a})^{\mathrm{T}}\right]$, where
$\mathbb{E}\left[\cdot\right]$ stands for the expectation operator and where
the upper script $(\cdot)^{\mathrm{T}}$ stands for the transpose operator
(later the adjoin operator for matrices). Hence, the distribution of the
initial condition is a Gaussian of mean $\chi^{a}$ and covariance matrix
$\mathbf{P}^{a}$, denoted by
$\mathcal{N}\left(\chi^{a},\mathbf{P}^{a}\right)$. The forecast state at a
time $T$, $\chi^{f}(T)=\mathbf{M}_{T\leftarrow 0}\chi^{a}$ provides an
approximation of the true state at time $T$. For linear dynamics and Gaussian
uncertainty, the forecast error $e^{f}(T)=\chi^{f}(T)-\chi^{t}(T)$ is a
Gaussian vector of zero mean and covariance matrix
$\mathbf{P}^{f}(T)=\mathbb{E}\left[e^{f}(e^{f})^{\mathrm{T}}\right](T)$, whose
dynamics is written as

$\partial_{t}e^{f}=\mathcal{M}(\partial e^{f}).$ (3)

The forecast-error covariance matrix is related to the analysis-error
covariance matrix by

$\mathbf{P}^{f}(T)=\mathbf{M}_{T\leftarrow
0}\mathbf{P}^{a}\left(\mathbf{M}_{T\leftarrow 0}\right)^{\mathrm{T}}.$ (4)

Equation (4) corresponds to the Kalman filter propagator of the error
covariance matrix, whose the particular dynamics is given by

$\frac{d\mathbf{P}^{f}}{dt}=\mathbf{M}\mathbf{P}^{f}+\mathbf{P}^{f}\mathbf{M}^{\mathrm{T}},$
(5)

integrated over the period $[0,T]$, starting from the initial condition
$\mathbf{P}^{f}(0)=\mathbf{P}^{a}$. Hence, the distribution of the forecast at
$T$ is a Gaussian of mean $\chi^{f}(T)$ and covariance matrix
$\mathbf{P}^{f}(T)$, that is
$\mathcal{N}\left(\chi^{f},\mathbf{P}^{f}(T)\right)$.

When observations, $y^{o}$, are available at $T$, so that
$y^{o}=\mathbf{H}\chi^{t}(T)+e^{o}$ where $\mathbf{H}$ denotes the linear
observation operator that maps a state (here the true state $\chi^{t}(T)$)
into the observation space, and $e^{o}$ denotes the observational error,
modeled as a Gaussian random vector of zero mean and covariance matrix
$\mathbf{R}$ and assumed uncorrelated from the forecast error $e^{f}(T)$ ;
then the Kalman analysis equations at time $T$ are written as

$\displaystyle\chi^{a}(T)$
$\displaystyle=\chi^{f}(T)+\mathbf{K}(y^{o}-\mathbf{H}\chi^{f}(T)),$ (6a)
$\displaystyle\mathbf{P}^{a}(T)$
$\displaystyle=(\mathbf{I}-\mathbf{K}\mathbf{H})\mathbf{P}^{f}(T),$ (6b)

where
$\mathbf{K}=\mathbf{P}^{f}(T)\mathbf{H}^{T}(\mathbf{H}\mathbf{P}^{f}(T)\mathbf{H}^{T}+\mathbf{R})^{-1}$
is the Kalman gain matrix. Eq. (6) characterizes the Gaussian distribution of
the analysis uncertainty, where Eq. (6a) is the equation of the mean, and Eq.
(6b) is the equation of the analysis-error covariance matrix at time $T$. Eq.
(4) together with Eq. (6) are the Kalman filter equations, that can be cycled
providing the analysis and forecast assimilation cycles, where Eq. (4) (Eq.
(6)) is the KF forecast (analysis) step.

Note that this contribution being focused on the uncertainties dynamics, with
the corresponding Kalman filter forecast of the covariance matrix Eq. (4), the
uncertainty related to the observations (present when considering cycled
assimilations) appears only through the statistics of the analysis error, and
more precisely through the covariance matrix of the analysis error
$\mathbf{P}^{a}=\mathbf{P}^{f}(0)$. Moreover, since we consider the model as
perfect, the uncertainty considered here only comes from the initial condition
and the boundary conditions.

Note also that the derivation of the Kalman filter forecast step Eq. (4)
follows the usual one encountered in the weather community, and can leads to
drawback when applied to LAMs without care e.g. a loss of variance at the
inlet of the domain. We detail this issue considering a transport in a 1D
limited-area domain in A. However, it results that when boundary conditions
are uncertain, it is necessary to reformulate the problem. That can be done in
the KF framework, by considering that the BCs are an uncertain command
following the terminology used in optimal control (see e.g. [42]).

###  2.2 Ensemble approximation of the Kalman filter forecast step

While the KF forecast step Eq. (4) is a simple algebraic formula, it fails to
apply in large systems because of its numerical cost: if $n$ denotes the
dimension of the vector representation of $\chi$, then the computational
complexity of Eq. (4) scales between $n^{2}$ and $n^{3}$ [44]. In terms of
integration cost, the KF requires $2n$ integrations of the model Eq. (1).

Hence, approximations for the KF are needed. For instance, in the Ensemble
Kalman filter (EnKF), the forecast error-covariance matrix is approximated by
its ensemble estimation.

$\widehat{\mathbf{P}^{f}}(t)=\frac{1}{N}\sum_{k}e_{k}e_{k}^{\mathrm{T}},$ (7)

with $e_{k}=\chi_{k}(t)-\widehat{\chi}(t)$ where
$\widehat{\chi}(t)=\frac{1}{N}\sum_{k}\chi_{k}(t)$ denotes the empirical mean
and $\left(\chi_{k}(t)\right)_{k\in[1,N]}$ is an ensemble of $N$ forecasts
[14]. This time, the numerical complexity scales with the number of ensemble
members $N$ and the size of the problem $n$: the numerical cost of an ensemble
of forecast is the cost of $N$ integrations of the model Eq. (1).

Note that the normalization by $N$ in Eq. (7) leads to a bias that decreases
as $1/N$. In EnKF framework, the normalization by $N-1$ is preferred, however
since we latter consider estimation from very large ensemble size, the
corrections of the estimators are not considered here, and we only consider
empirical mean estimations $\frac{1}{N}\sum_{k}(\cdots)$ as in Eq. (7).

In the ensemble Kalman filter forecast step applied for LAMs, the uncertainty
on the boundary condition can be introduced by adding perturbations at the
boundary for each member of the ensemble. This can be done from considering a
modeling of the multivariate covariance at the boundary or by considering
boundary condition derived from an ensemble on a larger domain [47].

While this contribution focuses on the forecast step, we can mention that if
ensemble computation of the KF analysis step Eq. (6) can be implemented
following different algorithm e.g. the EnKF version based on perturbation of
observation introduced by [7], the square root filter of [46] or the ETKF of
[2] and its variational implementation [17, 3] ; the KF forecast step is
mainly the same whatever the variant of the EnKF considered. Hence, for this
study dedicated to the exploration of the specification of the boundary
conditions for the PKF forecast step, the particular implementation the EnKF
has not impact and no assimilation experiment is conducted here.

The next section presents another approximation for the error-covariance
matrices.

###  2.3 Parametric formulation for the Kalman filter forecast step based on
VLATcov models

In the parametric approach, a covariance model is introduced,
$\mathbf{P}(\mathcal{P})$ where $\mathcal{P}$ denotes the set of parameters of
the covariance model, so to approximate the error covariance matrices. For
instance, the forecast-error covariance matrix $\mathbf{P}^{f}$, is
approximated as $\mathbf{P}(\mathcal{P}^{f})\approx\mathbf{P}^{f}$, where
$\mathcal{P}^{f}$ is a particular set of values for the parameters. The
parametric Kalman filter (PKF) dynamics aims to mimic the dynamics of Eq. (5)
relying on the dynamics of the parameters $\mathcal{P}^{f}$,

$\frac{d\mathcal{P}^{f}}{dt}=\mathcal{G}(\mathcal{P}^{f}),$ (8)

where $\mathcal{G}$ has to be determined from the particular dynamics of Eq.
(1), so that at any time $t$, $\mathbf{P}(\mathcal{P}^{f}(t))$ approximates
$\mathbf{P}^{f}(t)$ i.e.
$\mathbf{P}(\mathcal{P}^{f}(t))\approx\mathbf{P}^{f}(t)$. As for the EnKF, the
numerical complexity of the PKF prediction Eq. (8) scales as number of
parameters and the dimension of the problem $n$: the numerical cost of the PKF
represent the cost of few numerical integrations of the dynamics Eq. (1),
depending on the number of parameters needed for the covariance approximation.

Thereafter, since we deal with the forecast step of the PKF, the upper-script
f is dropped in the notation that concerns the forecast-error statistics.

This contribution will focus on the particular class of covariance model, so-
called VLATcov models, parameterized from two fields, defined below: the
variance field, $V$, and the local anisotropy tensor of the correlation
functions, $\mathbf{g}$ or $\mathbf{s}$. Hence, the set of parameters is given
by the couple $\mathcal{P}=(V,\mathbf{g})$ or $\mathcal{P}=(V,\mathbf{s})$, so
that a VLATcov model is written as $\mathbf{P}(V,\mathbf{g})$ or
$\mathbf{P}(V,\mathbf{s})$. For an error field $e$, the variance field is
defined as

$V=\mathbb{E}\left[e^{2}\right],$ (9)

When the error field is a differential random field, that is assumed from now,
the correlation function
$\rho(\mathbf{x},\mathbf{y})=\mathbb{E}\left[\varepsilon(\mathbf{x})\varepsilon(\mathbf{y})\right]$
is flat for $\mathbf{y}=\mathbf{x}$. Then, the local anisotropy at
$\mathbf{x}$ is defined as the local metric tensor $\mathbf{g}(\mathbf{x})$
(also denoted by $\mathbf{g}_{\mathbf{x}}$) which appears in the second-order
Taylor’s expansion

$\rho(\mathbf{x},\mathbf{x}+\delta\mathbf{x})\approx
1-\frac{1}{2}||\delta\mathbf{x}||_{\mathbf{g}_{\mathbf{x}}}^{2},$ (10)

where
$||\delta\mathbf{x}||_{\mathbf{g}_{\mathbf{x}}}^{2}=\delta\mathbf{x}^{\mathrm{T}}\mathbf{g}_{\mathbf{x}}\delta\mathbf{x}$
denotes the norm associated with the metric tensor $\mathbf{g}_{\mathbf{x}}$
that is the symmetric definite positive matrix
$[\mathbf{g}_{\mathbf{x}}]_{ij}=-\partial_{x^{i}x^{j}}^{2}\rho_{\mathbf{x}}$
where $\rho_{\mathbf{x}}(\delta\mathbf{x})$ stands for the local correlation
function. In a 1D domain of coordinate $x$, the metric tensor field is the
scalar field $\mathbf{g}=[g_{xx}]$.

In practice, the geometry of the local metric tensor is contravariant: the
direction of largest correlation anisotropy corresponds to the principal axes
of smallest eigenvalue for the metric tensor. Thus, it is useful to introduce
the local aspect tensor [11]

$\mathbf{s}(\mathbf{x})=\left(\mathbf{g}(\mathbf{x})\right)^{-1},$ (11)

where the superscript $(\cdot)^{-1}$ denotes the matrix inverse, and whose the
geometry goes as the correlation.

What makes the local metric tensor attractive is that this tensor is related
to the normalized error by (see e.g. [35])

$[\mathbf{g}_{\mathbf{x}}]_{ij}=\mathbb{E}\left[\partial_{x^{i}}\varepsilon\partial_{x^{j}}\varepsilon\right].$
(12)

Hence, the variance Eq. (9) and the anisotropy Eq. (12) can be computed from
an ensemble estimation: the variance field is estimated by

$\widehat{V}=\frac{1}{N}\sum_{k}\left(e_{k}(t)\right)^{2},$ (13)

with $e_{k}(t)=\chi_{k}(t)-\widehat{\chi}(t)$, from which derivatives of the
normalized error
$\varepsilon_{k}=\frac{1}{\sqrt{V}}\left(\chi_{k}(t)-\widehat{\chi}(t)\right)$
leads to the estimation of the upper triangular components of the metric

$\widehat{g}_{ij}=\frac{1}{N}\sum_{k}\partial_{x^{i}}\varepsilon_{k}\partial_{x^{j}}\varepsilon_{k},$
(14)

for $i\leq j$ (since $g_{ji}=g_{ij}$). While the PKF approach does not relies
on any ensembles, the ensemble estimations Eq. (13) and Eq. (14) can be used
to set the initial conditions for the parameters to ignite the assimilation
cycles, or to validate the PKF from the diagnosis of an EnKF.

An example VLATcov model is given by the heterogeneous Gaussian-like
covariance model [29]

$\mathbf{P}(V,\mathbf{s})(\mathbf{x},\mathbf{y})=\sqrt{V_{\mathbf{x}}V_{\mathbf{y}}}\frac{|\mathbf{s_{x}}|^{1/4}|\mathbf{s_{y}}|^{1/4}}{|\frac{1}{2}(\mathbf{s_{x}}+\mathbf{s_{y}})|^{1/2}}\exp\left(-\frac{1}{2}||\mathbf{x}-\mathbf{y}||^{2}_{[\frac{1}{2}(\mathbf{s_{x}}+\mathbf{s_{y}})]^{-1}}\right)$
(15)

where $|\cdot|$ denotes the matrix determinant.

When VLATcov models are used for the parametric approach, the dynamics of the
parameters Eq. (8) is deduced from the time derivative of Eq. (9) and Eq.
(12), and the dynamics of the error Eq. (3). Hence, in univariate statistics,
the PKF equations for VLATcov models consists in three equations. The first
equation corresponds to the prediction of the mean, that is Eq. (1) in the
linear case. The second equation is the dynamics of the variance, $V$, which
is deduced from

$\partial_{t}V=2\mathbb{E}\left[e\partial_{t}e\right],$ (16a) where replacing
the trend of the error Eq. (3), will leads to the dynamics of $V$
$\partial_{t}V=2\mathbb{E}\left[e\mathcal{M}(\partial e)\right].$ (16b)

The third equation is the dynamics of the anisotropy, $g$, whose components
evolve as

$\partial_{t}g_{ij}=\mathbb{E}\left[\partial_{t}\left(\partial_{x^{i}}\varepsilon\partial_{x^{j}}\varepsilon\right)\right],$
(17)

and from which we can deduce the dynamics of the aspect tensor
$\mathbf{s}=\mathbf{g}^{-1}$ with
$\partial_{t}\mathbf{s}=-\mathbf{s}(\partial_{t}\mathbf{g})\mathbf{s}$. The
dynamics of the variance Eq. (16b) and of the anisotropy Eq. (17) can be
simplified e.g. by considering the commutation between the expectation and
partial derivatives [36].

Note that for non-linear dynamics, the PKF forecast step for VLATcov models
stands as (see [36] for details)

$\displaystyle\partial_{t}\mathbb{E}\left[\chi\right]$
$\displaystyle=\mathcal{M}(\partial\mathbb{E}\left[\chi\right])+\mathcal{M}^{\prime\prime}(\partial\mathbb{E}\left[\chi\right])\left(\mathbb{E}\left[\partial
e\otimes\partial e\right]\right),$ (18a) $\displaystyle\partial_{t}V$
$\displaystyle=2\mathbb{E}\left[e\partial_{t}e\right],$ (18b)
$\displaystyle\partial_{t}g_{ij}$
$\displaystyle=\mathbb{E}\left[\partial_{t}\left(\partial_{x^{i}}\varepsilon\partial_{x^{j}}\varepsilon\right)\right],$
(18c)

where $\mathcal{M}^{\prime}$ and $\mathcal{M}^{\prime\prime}$ with two linear
operators : the former (the latter) refers to the tangent-linear model (the
Hessian), and both are computed with respect to the mean state
$\mathbb{E}\left[\chi\right]$ ; $\partial e\otimes\partial e$ denotes the
tensor product of the partial derivatives with respect to the spatial
coordinates, i.e., terms such as $\partial^{k}e\partial^{m}e$ for any positive
integers $(k,m)$ ; and where the trends than appears in the right-hand side of
Eq. (18) should be formulated from the dynamics e.g. using
$\partial_{t}e=\mathcal{M}^{\prime}(\partial\mathbb{E}\left[\chi\right])(\partial
e)$ for the trend of $e$.

The computation of dynamical equations Eq. (18) for the mean, the variance $V$
and the anisotropy $\mathbf{g}$ (or $\mathbf{s}$) can be performed using a
computed algebra system. To do so, the open source Python toolbox SymPKF has
been introduced [30, 36], which computes the dynamics of the parameters and
renders a numerical code to facilitate the numerical exploration of the PKF
approach. Another way to simplify the computation of the parameters dynamics
is to identify the contribution of each physical process present in Eq. (1)
following a splitting strategy [32, 36]. Thereafter, the dynamics of the
VLATcov parameters is computed by using SymPKF and the interested reader is
referred to the Jupyter notebooks that are provided as a supplementary
material to this contribution [41].

Compared to the ensemble implementation of the KF (e.g. the EnKF or the ETKF),
the contribution of the PKF is that it replaces the computational cost for the
ensemble prediction step by the numerical cost of the time integration of the
system Eq. (18). More precisely, in univariate statistics, the PKF based on
the VLATcov model scales as the number of independent components in
$\mathbf{g}$ (the number of coefficients in the upper triangle) plus one for
the variance field: in a univariate over a 1D (3D) domain, this represents 2
(7) times the cost of one model forecast (which scales itself with the
dimension $n$). This numerical cost of the PKF is then competitive to the
ensemble method where it often needs dozen of members. However in multivariate
statistics, the numerical cost of the PKF scales as the square of the number
of fields, which is a strong limitation e.g. a multivariate assimilation in
air quality should consider hundreds of chemical species [39], but in
practice, only a few species are assimilated, which makes the PKF interesting
for these applications (in CAMS regional air quality production 2.40 [8], the
univariate 3DVar system of MOCAGE is used for the separated assimilation of
ozone, nitrogen dioxide, sulphur dioxide, and fine particulate matter PM2.5
and PM10, following a configuration similar to the one used for MACII detailed
by [23]).

While it is not the purpose of the present work, we can mention that the
analysis step in the PKF consist in a sequential update of the parameters
during the assimilation of observations i.e. the update of the variance and of
the anisotropy for the VLATcov model after the assimilation of each
observation. For instance, at the leading order, the update of the aspect
tensor field due to the assimilation of an in situ observation is written as

$\mathbf{s}^{a}=\frac{V^{a}}{V^{f}}\mathbf{s}^{f},$ (19)

where $V^{a}$ is the analysis-error variance field that results from the
assimilation of the observation (see update equations Eq.(3) for the variance
and Eq.(12) for the aspect tensor in [35] for details). In the PKF, a
sequential processing of batches of observation can be consider to parallelize
the assimilation, similarly to the procedure employed in the EnKF analysis
step [18].

The PKF based on the VLATcov model is illustrated in the next sections for two
dynamics which give an explicit form for $\mathcal{M}$ in Eq. (1).

###  2.4 Illustration of the PKF forecast step for simple dynamics

The transport and the diffusion equations are considered so to detail the
dynamics of the variance and the anisotropy for the PKF applied for VLATcov
models. Both dynamics play over a 1D periodical domain of coordinate $x$, so
that the dynamics is an evolution equation without boundary conditions. The
PKF for the transport has been already detailed in previous contributions (see
e.g. [36] ), and is recapped here for self consistency. The PKF for the
heterogeneous diffusion equation is original, and extends the derivation of
the PKF applied to the homogeneous diffusion equation [32].

####  2.4.1 PKF prediction applied on a transport equation

The transport equation of a scalar field $c(t,x)$ by a stationary velocity
field $u(x)$ is written as

$\partial_{t}c+u\partial_{x}c=0.$ (20)

In this example, and by identification with Eq. (1), $c$ stands for $\chi$
while $\mathcal{M}(c,\partial c)=-u\partial_{x}c$. This kind of equation
appears for instance in the prediction of the concentration of a chemical
specie as in chemical transport models.

The computation of the PKF dynamics for Eq. (20) using SymPKF leads to the
system

$\displaystyle\partial_{t}c$ $\displaystyle=-u\partial_{x}c,$ (21a)
$\displaystyle\partial_{t}V_{c}$ $\displaystyle=-u\partial_{x}V_{c},$ (21b)
$\displaystyle\partial_{t}s_{c,xx}$
$\displaystyle=-u\partial_{x}s_{c,xx}+2s_{c,xx}\partial_{x}u,$ (21c)

where the anisotropy is represented by the aspect tensor $\mathbf{s}=s_{c,xx}$
in 1D domain. The PKF dynamics Eq. (21) is a system of three uncoupled partial
derivative equation similar to the one first found by [9]. This system
represents the dynamics of the mean state $\mathbb{E}\left[c\right]$, Eq.
(21a), where the expectation operator has been removed for the sake of
simplicity ; the transport of the variance, Eq. (21b) ; and the transport of
the anisotropy Eq. (21c), where an additional a source term of anisotropy
appears, that is due to the shear by the flow. Compared with an ensemble
approach, the PKF approach allows for an understanding of the dynamics and the
physics of the uncertainty.

This example shows that the PKF does not need solving the full state dimension
covariance matrix, but only the dynamics of the parameters used in the
modeling of the covariance matrix: the variance and the anisotropy in the
VLATcov model considered here which evolves from Eq. (21b) and Eq. (21c)
respectively. However, the parallelization of the PKF is less obvious than for
the parallel computation of an ensemble of forecast, since it requires to
parallelize the coupled system Eq. (21).

Note that the lower script notation c for $V_{c}$ and c,xx for $s_{c,xx}$
corresponds to the notation automatically rendered by SymPKF when processing
the dynamics Eq. (20) at a symbolic level. This labelling for the parameters
has been introduced when multiple fields are present e.g. in multivariate
dynamics. While this contribution only address univariate dynamics, the
notation is kept here so to facilitate the comparison with the output of
SymPKF and also because another important dynamics is discussed: the diffusion
equation, which is now presented.

####  2.4.2 PKF prediction applied on a diffusion equation

The diffusion equation of a scalar field $f(t,x)$ and of diffusion coefficient
$D(x)$,

$\partial_{t}f=\partial_{x}\left(D\partial_{x}f\right),$ (22)

is now considered. This kind of equation appears for instance in the
prediction of electron density $f$ of the Earth radiation belts and results
from a Hamiltonian formalism applied on a typical Boltzmann equation, where a
Fokker-Planck operator is introduced to evaluate physical interactions
responsible for changing particles trapping state [10]. In the radiation
belts, the typical spatial coordinates system $x$ in Eq. (22) stands in this
case for a combined spatial and physical quantities e.g. the energy of the
electrons. The diffusion equation is also important in the modeling of
atmospheric boundary layer where it represents the effect of the turbulence
[45]. In this example, and by identification with Eq. (1), $f$ stands for
$\chi$ while $\mathcal{M}(f,\partial
f)=\partial_{x}\left(D\partial_{x}f\right)$.

The computation of the PKF dynamics for Eq. (22) can be performed using
SymPKF. However, because of the second order derivative, the dynamical system
makes appear an unknown term
$\mathbb{E}\left[\varepsilon_{f}\partial^{4}_{x}\varepsilon_{f}\right]$, not
determined from $f$, $V_{f}$ and $s_{f,xx}$ (see B). An analytical closure has
been proposed for 1D domains which states as [32]

$\mathbb{E}\left[\varepsilon_{f}\partial^{4}_{x}\varepsilon_{f}\right]=3g_{f,xx}^{2}-2\partial^{2}_{x}g_{f,xx}$
(23a) when written in metric tensor or
$\mathbb{E}\left[\varepsilon_{f}\partial^{4}_{x}\varepsilon_{f}\right]=\frac{2\partial^{2}_{x}s_{f,xx}}{s_{f,xx}^{2}}+\frac{3}{s_{f,xx}^{2}}-\frac{4\left(\partial_{x}s_{f,xx}\right)^{2}}{s_{f,xx}^{3}}$
(23b) in aspect tensor, which leads to the PKF dynamics
$\displaystyle\partial_{t}f$
$\displaystyle=D\partial^{2}_{x}f+\partial_{x}D\partial_{x}f,$ (24a)
$\displaystyle\partial_{t}V_{f}$
$\displaystyle=-\frac{2DV_{f}}{s_{f,xx}}+D\partial^{2}_{x}V_{f}-\frac{D\left(\partial_{x}V_{f}\right)^{2}}{2V_{f}}+\partial_{x}D\partial_{x}V_{f},$
(24b) $\partial_{t}s_{f,xx}=D\partial^{2}_{x}s_{f,xx}+4D\\\
-\frac{2D\left(\partial_{x}s_{f,xx}\right)^{2}}{s_{f,xx}}-\frac{2Ds_{f,xx}\partial^{2}_{x}V_{f}}{V_{f}}+\\\
\frac{D\partial_{x}V_{f}\partial_{x}s_{f,xx}}{V_{f}}+\frac{2Ds_{f,xx}\left(\partial_{x}V_{f}\right)^{2}}{V_{f}^{2}}-2s_{f,xx}\partial_{x}^{2}D+\\\
2\partial_{x}D\partial_{x}s_{f,xx}-\frac{2s_{f,xx}\partial_{x}D\partial_{x}V_{f}}{V_{f}},$
(24c)

where in this dynamical systems, the expected value $\mathbb{E}\left[f\right]$
in Eq. (24a) is replaced by $f$ for the sake of simplicity. The dynamics Eq.
(24) makes appear the effect of the transport due the heterogeneity of the
diffusion coefficient which implies a flow of velocity $-\partial_{x}D$, and
leads to the same PKF transport dynamics Eq. (21) as discussed for Eq. (20) in
the particular case where $u=-\partial_{x}D$. The other terms in Eq. (24) are
related to the second-order derivative term $D\partial_{x}^{2}f$, which
couples the dynamics of the variance and of the anisotropy.

In term of metric, the closed Eq. (24) reads as

$\displaystyle\partial_{t}f$
$\displaystyle=D\partial^{2}_{x}f+\partial_{x}D\partial_{x}f,$ (25a)
$\displaystyle\partial_{t}V_{f}$
$\displaystyle=-2DV_{f}g_{f,xx}+D\partial^{2}_{x}V_{f}-\frac{D\left(\partial_{x}V_{f}\right)^{2}}{2V_{f}}+\partial_{x}D\partial_{x}V_{f},$
(25b) $\partial_{t}g_{f,xx}=-4Dg_{f,xx}^{2}+D\partial^{2}_{x}g_{f,xx}+\\\
\frac{2Dg_{f,xx}\partial^{2}_{x}V_{f}}{V_{f}}+\frac{D\partial_{x}V_{f}\partial_{x}g_{f,xx}}{V_{f}}-\frac{2Dg_{f,xx}\left(\partial_{x}V_{f}\right)^{2}}{V_{f}^{2}}+\\\
2g_{f,xx}\partial_{x}^{2}D+2\partial_{x}D\partial_{x}g_{f,xx}+\frac{2g_{f,xx}\partial_{x}D\partial_{x}V_{f}}{V_{f}}.$
(25c)

Compared with the use of an ensemble estimation for the error statistics, the
dynamics Eq. (25) appears rather complicated, especially the dynamics of the
metric Eq. (25c). However, this coupled system is simple to solve numerically
e.g. considering a finite difference scheme and automatic code generation.
Moreover, the numerical cost of Eq. (25) is of the order of three times the
cost the original heterogeneous diffusion equation Eq. (22) while it provides
an accurate estimation of the error statistics (see numerical experiment in
Section 5.3), where a large ensemble is needed to reach this accuracy and
avoid the sampling noise whose magnitude, for an ensemble of$N_{e}$ members,
scales in $1/\sqrt{N_{e}}$ from the central limit theorem (CLT).

Until now, PKF dynamics for the heterogeneous diffusion equation has been
evaluated on periodic domain only, while bounded domains are often needed,
e.g. in radiation belts predictions where the energy of electrons are limited,
or in atmospheric boundary layer where the ground is a limit of the domain.
The next section addresses how to specify the boundary conditions for the PKF
dynamics.

##  3 Specification of the boundary conditions for the PKF forecast step

This section tackles the specification of the boundary conditions for the PKF
by considering two usual kind of conditions: the Dirichlet and the Neumann
conditions. We consider the particular case of the semi-bounded 1D domain
$[0,\infty)$, and focus on the boundary $x=0$. Then we extend to boundary
conditions of an arbitrary domain $\Omega$ of frontier $\partial\Omega$.

###  3.1 Dirichlet BCs

A Dirichlet condition at the boundary consists in specifying the value of the
fields at $x=0$, that is $\chi(t,x=0)=\chi_{0}(t)$.

This conditions is used for the dynamics of the mean in the PKF, but it
remains to specify the boundary conditions for the variance and the
anisotropy.

Therefore the Dirichlet condition implies that the error field must also
verifies a Dirichlet condition i.e. $e(t,x=0)=e_{0}(t)$. The expectation of
the error field at $x=0$ is zero by definition, and of variance
$V_{0}(t)=\mathbb{E}\left[e_{0}(t)^{2}\right]$. Hence, the variance field must
also verify a Dirichlet condition i.e. $V(t,x=0)=V_{0}(t)$.

So for a 1D bounded domain, the Dirichlet condition on the dynamics implies to
specify a Dirichlet condition on the variance and on the anisotropy. This
result extends for an arbitrary domain $\Omega$ where this time, the boundary
conditions for the variance and the anisotropy are Dirichlet conditions on the
frontier $\partial\Omega$.

In case where the bounded domain is nested within a larger domain where
uncertainty is known from a PKF dynamics, then the variance and the anisotropy
at the boundary can be set from the variance and the anisotropy known in the
larger domain. When the uncertainty at large scale is featured from an
ensemble of forecasts, the statistics at the boundary should be set as the
statistics estimated from the ensemble of large scale forecasts at the
boundary points e.g. for VLATcov models, the variance and the anisotropy can
be estimated from the ensemble of large scale forecasts from Eq. (13) and Eq.
(14) respectively.

Hence, Dirichlet condition in case of nested models easily extends in 2D and
3D domains where it remains to specify the variance and the anisotropy of the
local area model from the variance and the anisotropy of the coupling model.

###  3.2 Neumann BCs

Neumann conditions at the boundaries are written as null fluxes i.e.
$\partial_{x}\chi(t,x=0)=0$. This implies that the error field must also
verifies a Neumann condition i.e. $\partial_{x}e(t,x=0)=0$. Again, we are
looking for the boundary conditions for the variance and the anisotropy.

The condition on the variance is deduced from the Taylor expansion of the
error at the vicinity of $x=0$ as follows. The expectation of the square of
the second order expansion of the error

$e(t,\delta x)=e(t,0)+\frac{1}{2}\partial_{x}^{2}e(t,0)\delta
x^{2}+\mathcal{O}(\delta x^{3}),$

leads to the local expansion of the variance
$V(t,x)=\mathbb{E}\left[e^{2}\right](t,x)$,

$V(t,\delta x)=V(t,0)+\mathbb{E}\left[e\partial_{x}^{2}e\right](t,0)\delta
x^{2}+\mathcal{O}\left(\delta x^{4}\right).$

As the local Taylor expansion of the variance field at $x=0$, this implies
that the first order derivative is null, i.e.

$\partial_{x}V(t,0)=0,$ (26a) which means that the condition in variance at
the boundary $x=0$ follows a Neumann condition.

For the anisotropy, the Neumann condition on the variance,Eq. (26a), implies
that the metric tensor
$g(t,x)=\mathbb{E}\left[\left(\partial_{x}\varepsilon\right)^{2}\right](t,x)$
simplifies as
$g(t,0)=\frac{1}{V(t,0)}\mathbb{E}\left[\left(\partial_{x}e(t,0)\right)^{2}\right]$.
Then the Neumann condition on $e$, $\partial_{x}e(t,0)=0$ i.e. implies that
the condition for the metric is a Dirichlet condition,

$g(t,x=0)=0.$ (26b)

Note that the Dirichlet condition for the metric Eq. (26b) is equivalent to an
infinite in aspect tensor at the boundary i.e. $s(t,0)=+\infty$. Therefore,
the PKF formulated in aspect tensor is singular near the boundaries where the
error field verifies a Neumann condition. Hence, for a Neumann condition on
the error, it is preferable to consider the PKF formulation in variance/metric
tensor rather than in variance/aspect tensor (while there is no preference on
unbounded domains e.g. in periodic domains or flat spaces).

Hence, the Neumann condition on the error for a 1D domain leads, for the PKF
dynamics, to a Neumann condition in variance and a Dirichlet condition in
metric. This extends to a 2D or 3D domain $\Omega$ where this time the Neumann
condition for the variance states as a null flux along the normal direction of
the frontier $\partial\Omega$ of the domain. The Dirichlet condition for the
metric reads equivalently $g(t,\mathbf{x})=0$ for
$\mathbf{x}\in\partial\Omega$, but a weaker condition could be introduced
where the tangential components of the metric at the boundary are not zero
(not addressed here).

Now that the boundary conditions for the PKF have been theoretically specified
for the Dirichlet and the Neumann conditions, a numerical validation as well
as a comparison with the usual EnKF approach is introduced. But to do so, it
is necessary to specify an appropriate setting for the boundary of the EnKF,
as discussed in the next section.

##  4 Methodology for validating the PKF from an EnKF simulations, and
Specification of the BCs for EnKF

###  4.1 Methodology for the validation of the PKF from an ensemble method

In the previous sections we specified the PKF equations as well as the
boundary conditions it needs depending on the boundary condition for the error
at the boundary. Note that for linear dynamics as for the transport and the
diffusion that are considered here, the dynamics of the error statistics is
independent from the dynamics of the mean, so hereafter we will not be
considering the dynamics of the mean but rather focusing on the error
statistics. For instance, the PKF dynamics for the error statistics for the
transport with a positive wind on a bounded 1D domain and Dirichlet condition
in $x=0$, Eq. (21) is written as

$\displaystyle\partial_{t}V_{c}$ $\displaystyle=-u\partial_{x}V_{c},$ (27a)
$\displaystyle V_{c}(0,x)$ $\displaystyle=V_{c}^{0}(x),$ (27b) $\displaystyle
V_{c}(t,0)$ $\displaystyle=b_{V_{c}}(t),$ (27c)
$\displaystyle\partial_{t}s_{c,xx}$
$\displaystyle=-u\partial_{x}s_{c,xx}+2s_{c,xx}\partial_{x}u,$ (27d)
$\displaystyle s_{c,xx}(0,x)$ $\displaystyle=s_{c,xx}^{0}(x),$ (27e)
$\displaystyle s_{c,xx}(t,0)$ $\displaystyle=b_{s_{c,xx}}(t),$ (27f)

where $b_{V_{c}}$ and $b_{s_{c,xx}}$ stands respectively for the error
variance and anisotropy time functions prescribed at the inlet of the domain
in $x=0$, and where $V_{c}^{0}$ and $s_{c,xx}^{0}$ are the initial condition
fields.

Now we are interested to verify that the resulting PKF dynamics, such as Eq.
(27), is really able to reproduce the dynamics of the error statistics as an
ensemble method would provide. But then, it is necessary to create an ensemble
of forecast with an appropriate specification of the boundary error. Again,
for the transport with Dirichlet condition, it remains to perform an ensemble
of integration of the dynamics

$\displaystyle\partial_{t}e_{c}$ $\displaystyle=-u\partial_{x}e_{c},$ (28a)
$\displaystyle e_{c}(0,x)$ $\displaystyle=e_{c}^{0}(x),$ (28b) $\displaystyle
e_{c}(t,0)$ $\displaystyle=b_{e_{c}}(t),$ (28c)

where $b_{e_{c}}$ and $e_{c}^{0}$ stands respectively for the error time
function prescribed at the inlet of the domain in $x=0$, and the initial error
field.

In real application where the inflow at the inlet is related to physical
structures e.g. a mid-latitude vortex entering in a LAM, the error is
continuous and differentiable in time and space meaning this imposes a
constraint of coherence between the error field inside the domain and its
boundary condition at the inlet. Then, the resulting error statistics should
also be continuous and differentiable in time and space: this imposes a
constraint of coherence between the variance (or the anisotropy) field inside
the domain and its boundary condition at the inlet.

Note that when a covariance model is used in EnKF to perturb boundary
conditions, an assumption of stationarity of the error statistics is often
encountered, e.g. by considering a constant time correlation in an auto-
regressive time process [47] (similar to the one introduced in Appendix A).
However for physical phenomenons entering in a LAM, the correlation time scale
is expected to evolve accordingly to the spatial anisotropy, because the
spatial anisotropy is governed by the dynamics and the observational network
[9, 5, 35]. It would therefore be interesting to illustrate how the PKF
performs for heterogeneous and non-stationary error statistics at the
frontier.

Thus, we consider the following validation framework where the error
statistics are given by a prescribed initial condition together with its
boundary conditions which leads to smooth evolution of the error statistics,
without any constraint of homogeneity and stationarity of the statistics. For
instance, when considering the transport, this means that we consider the
situation where the error statistics are given by a prescribed initial
condition $\left(V_{c}^{0}(x),s_{c,xx}^{0}(x)\right)$ together with its
boundary conditions $\left(b_{V_{c}}(t),b_{s_{c,xx}}(t)\right)$ which leads to
smooth fields $\left(V_{c}(t,x),s_{c,xx}(t,x)\right)$. But, for the ensemble
use for the validation, this implies to populate an ensemble of $N_{e}$
initial conditions $\left({e_{c}}_{k}(0,x)\right)_{k\in[1,Ne]}$ together with
their boundary conditions $\left({b_{e_{c}}}_{k}(t)\right)_{k\in[1,Ne]}$ such
that at a given time $t$ and position $x$, the resulting error variance
(anisotropy) is $V_{c}(t,x)$ ($s_{c,xx}(t,x)$). The problem that arises is how
to specify the auto-correlation time scale of the boundary perturbation in
order to obtain the desired error statistics
$\left(V_{c}(t,x),s_{c,xx}(t,x)\right)$.

Note that ensemble forecasting for the diffusion equation with Neumann
boundary conditions corresponds to an initial value problem where each member
is integrated from an initial condition that verifies the Neumann conditions,
and without additional perturbation of the boundary along the time integration
(the metric at the boundary is constant equal to zero).

It results that, the specification of the auto-correlation time scale only
concerns the Dirichlet conditions. In what follows, the specification of the
time auto-correlation is first presented for an arbitrary evolution equation,
then it is applied for the transport and for the diffusion equation.

###  4.2 Specification of the auto-correlation time-scale of the BC
perturbations for ensemble of forecast

The problem faced here is that, with boundary conditions being time-series, it
has an associated time-scale. However, the metric tensor $g_{xx}$ is related
to the spatial length-scale of the perturbation as denoted by the index xx. In
order to specify the boundary condition of the metric tensor field, we need to
find an equation linking, on the boundaries, the spatial metric tensor
$g_{xx}$ with the time-scale used to generate the perturbation.

Similarly to the spatial metric tensor Eq. (12), the temporal metric tensor
$g_{tt}$ that characterize the auto-correlation of a smooth centered random
field, $\eta(t)$, depending on the time, and of variance
$V_{\eta}(t)=\mathbb{E}\left[\eta(t)^{2}\right]$, is defined by

$\mathbf{g}_{tt}(t)=\mathbb{E}\left[\partial_{t}\left(\frac{\eta(t)}{V_{\eta}(t)}\right)\partial_{t}\left(\frac{\eta(t)}{V_{\eta}(t)}\right)\right].$
(29)

This temporal metric tensor is directly related with the time-scale of the
perturbation. In 1D $g_{tt}=\frac{1}{L_{t}^{2}}$ with $L_{t}$ the auto-
correlation time-scale.

Without loss of generality, the boundary $x=0$ is considered, and the goal is
to characterize the temporal metric tensor $g_{tt}(t,x=0)$. If $\eta(t)$
denotes the random error at $x=0$, then by continuity, the error and the
random forcing verify $e(t,x=0)=\eta(t)$. Then, it results that the variances
verify $V_{\eta}(t)=V(t,x=0)$, and the temporal metric tensor reads as

$\mathbf{g}_{tt,x=0}(t)=\mathbb{E}\left[\partial_{t}\varepsilon(t,\mathbf{x}=0)\partial_{t}\varepsilon(t,\mathbf{x}=0)\right],$
(30a) where $\varepsilon=e/\sqrt{V}$ is the normalized error associated with
the spatial error $e$. While Eq. (30a) only holds at the boundary $x=0$, the
spatio-temporal smoothness of $e$ implies a link between the temporal metric
at the boundary and the spatial metric within the domain, which results from
the dynamics of the error Eq. (3) at $x=0$:
$\partial_{t}e(t,x=0)=\mathcal{M}(e,\partial e)(t,x=0)$. In particular, the
temporal metric reads as (see C)
$g_{tt}\underset{x=0}{=}\frac{1}{V}\mathbb{E}\left[\left(\mathcal{M}(\varepsilon\sqrt{V},\partial(\varepsilon\sqrt{V}))\right)^{2}\right]-\frac{1}{4V^{2}}\left(\partial_{t}V\right)^{2},$
(30b) where the terms
$\mathbb{E}\left[\left(\mathcal{M}(\varepsilon\sqrt{V},\partial(\varepsilon\sqrt{V}))\right)^{2}\right]$
and $\partial_{t}V$ can make appear the spatial metric field $g_{xx}$ at
$x=0$.

One pitfall is that equation Eq. (30b) may be complicated, and can contain
unknown terms such as
$\mathbb{E}\left[\varepsilon_{f}\partial^{4}_{x}\varepsilon_{f}\right]$
encountered for the heterogeneous diffusion dynamics in section 2.4.2. The
next two sub-sections will detail the link between the temporal and the
spatial metrics for the transport and for the diffusion.

###  4.3 Dirichlet BC for ensemble forecasting of the positive velocity
transport equation

To illustrate the relation between the temporal and the spatial metric tensor,
the transport equation Eq. (20) is now considered.

Following the theoretical derivation of the temporal metric Eq. (30b), SymPKF
is used to derive the relation between the temporal anisotropy tensor of the
error and the spatial anisotropy tensor of the error which yields

$g_{c,tt}\underset{x=0}{=}u^{2}g_{c,xx}+\frac{u^{2}\left(\partial_{x}V_{c}\right)^{2}}{4V_{c}^{2}}+\frac{u\partial_{t}V_{c}\partial_{x}V_{c}}{2V_{c}^{2}}+\frac{\left(\partial_{t}V_{c}\right)^{2}}{4V_{c}^{2}}.$
(31)

This spatio-temporal consistency for the temporal and spatial statistics is
difficult to interpret physically without approximations. However, under the
assumptions of local homogeneity ($\partial_{x}V_{c}=0$) and of stationarity
for the variance $\partial_{t}V_{c}=0)$, Eq. (31) reads as

$g_{c,tt}\underset{x=0}{=}u^{2}g_{c,xx},$ (32)

which is physically interpretable since Eq. (32), written in time-scale and
length-scale, reads as $L_{t}=\frac{L_{x}}{u}$: the usual rule relating time
and space in a transport. Later, the numerical investigation will consider Eq.
(32) as an approximation of the true time-scale even when assumptions leading
to Eq. (32) are not verified.

Note that Eq. (32) can be obtained when considering that the dynamics of the
variance Eq. (21b) applies at the boundary, leading to replace the trend of
the variance by $\partial_{t}V_{c}=-u\partial_{x}V_{c}$ in Eq. (31) so to
obtain Eq. (32).

To conclude this paragraph, the ensemble forecasting under Dirichlet boundary
conditions and applied to the transport equation, remains to populate an
ensemble of boundary perturbations with a prescribed temporal variance and an
auto-correlation time scale given by Eq. (32).

We proceed in the same way for the diffusion equation.

###  4.4 Dirichlet BCs for ensemble forecasting of the diffusion equation

To continue going towards more and more realistic modeling, the heterogeneous
diffusion equation Eq. (22) is now considered to compute the spatio-temporal
link Eq. (30b) in the diffusion case.

From a derivation detailed in D, the auto-correlation time scale of boundary
perturbation can be related to the spatial error correlation length-scale by
the proxy

$g_{f,tt}(t,x)\approx 3D(x)^{2}g_{f,xx}(t,x).$ (33)

Note that Eq. (33) is an equality when the variance and the diffusion fields
are homogeneous, and when the variance is stationary at the boundary.

Hence, as for the transport, the ensemble forecasting under Dirichlet boundary
conditions, and applied to the diffusion equation, remains to populate an
ensemble of boundary perturbations with a prescribed temporal variance and an
auto-correlation time scale given by Eq. (33).

We are now ready to validate the PKF approach from an ensemble validation
designed to produce desired error statistics.

##  5 Numerical investigation

The goal of the numerical investigation is to validate the PKF on a bounded
domain as well as the equations developed in Section 4, by comparing the PKF
dynamics with an ensemble simulation.

###  5.1 Setting for the numerical experiments

For this investigation three different settings are considered. All
experiments take place on a 1D bounded domain $x\in[0,\Lambda]$. For the first
one, the transport equation Eq. (20) is considered with Dirichlet boundary
condition at $x=0$ and free boundary at $x=\Lambda$. For the second setting,
the heterogeneous diffusion equation Eq. (22) is considered with Dirichlet
boundary conditions at both boundaries $x=0$ and $x=\Lambda$. For the third
setting, the same diffusion equation is considered but this time with Neumann
boundary conditions at $x=0$ and $x=\Lambda$.

The transport and the diffusion being linear, the dynamics of the mean is the
same for the PKF and for the EnKF. Hence, without loss of generality, to focus
on the validation of the error statistics, the mean state is not considered in
the following (the reader can consider the mean state as constant). Then, the
ensemble of forecast is equivalent to the forecasts of an ensemble of
perturbations $(e_{k})_{k\in[1,N_{e}]}$, with appropriate boundary conditions.

Each time the variance, Eq. (21b) and Eq. (24b), and anisotropy tensor, Eq.
(21c) and Eq. (24c), produced by the PKF dynamics are compared with the
variance and anisotropy tensor diagnosed from an ensemble of $N_{e}=6400$
forecasts. The size of the ensemble considered here, $6400$, is quite unusual
compared with operational ensembles that often count dozens of members ; but
the reason is that we want to avoid sampling noise to be able to verify that
the PKF dynamics of the error statistics reproduces the EnKF considered as the
reference. Hence, if a difference appears between the EnKF and the PKF, then
it is less likely to be related to sampling noise except when the magnitude of
the difference falls under the confident interval deduced from the CLT.

The domain is discretized in $n=241$ grid points and the spatial derivative
operator $\partial_{x}$ is discretized with a centered finite difference
scheme leading to a second order of consistency. The temporal discretization
scheme varies with each experiment and is detailed in each sections.

###  5.2 Application to the transport equation

Figure 1: Heterogeneous velocity field considered for the numerical simulation
of the transport dynamics. Figure 2: Sample of a generated perturbation split
into an initial condition and a boundary condition that are smoothly
connected. Figure 3: Comparison of the forecast-error variance (left column)
and normalized length-scale (right column) fields dynamics for the
heterogeneous advection equation on a 1D bounded domain with Dirichlet
boundary conditions at $x=0$ and open boundary condition at $x=\Lambda$. The
results are shown for times $t=0$, $t=0.25T_{\mathrm{adv}}$ and
$t=1.5T_{\mathrm{adv}}$.  Figure 4: Time evolution of the forecast-error
variance (a) and normalized length-scale (b) at $x=0.5\Lambda$, for the
advection equation with Dirichlet boundary conditions.

In this experiment setting, the transport equation Eq. (20) is considered. The
velocity wind for the simulation is set as the heterogeneous stationary field
$u(x)=1+\frac{1}{4}\sin(\frac{2\pi}{\Lambda}x)$ shown in Fig. 1.

The temporal discretization scheme used for the ensemble simulation as well as
the PKF dynamics is a Runge-Kutta scheme of order 4 with a fixed time-step
$dt\approx 4.10^{-3}$. The simulation is conducted from time $t=0$ until
$t_{end}=2T_{adv}$ with the advection time scale
$T_{adv}=\frac{\Lambda}{u_{max}}$ .

In order to generate a coherent set of perturbations for the ensemble
simulation i.e. an initial condition and a boundary condition that are
smoothly connected, an extended domain $[-u(0)t_{end},\Lambda]$ is created
from the union of the physical domain $[0,\Lambda]$ and the time window
$[0,t_{end}]$ brought back to a virtual physical extension of the domain by
multiplying with $-u(0)$.

Then on this extended domain a variance field, $V_{0}$ and a length-scale
field $L_{0}$ are defined which will be used to generate the perturbations.
For this experiment the fields $V_{0}$ and $L_{0}$, that constitute the PKF
initial and boundary conditions, are set as follows. The initial variance is
set homogeneous and equal to 1 over the physical domain $V_{0}(t=0,x)=1$ and
the boundary variance is set to the periodical function
$V_{0}(t,x=0)=\frac{5}{4}-\frac{1}{4}cos(\frac{2\pi}{T_{adv}}t)$. Like for the
variance, the initial length-scale is set homogeneous and equal to 10% of the
domain length $L_{0}(t=0,x)=L_{h}=0.1\Lambda$ and the boundary length-scale is
set to the periodical function
$L_{0}(t,x=0)=0.1\Lambda(\frac{3}{4}+\frac{1}{4}cos(\frac{2\pi}{T_{adv}}t))$.

This setting for the variance and the length-scale is chosen so to represent a
typical behaviour encountered in numerical weather forecasting, where large
scale are more predictable than small scales, which is also the case in
radiation belts dynamics forecasting.

Using Eq. (15) and the relation between the length-scale and the anisotropy
tensor in 1D, $s_{0}=L_{0}^{2}$, the covariance matrix,
$\mathbf{P}_{0}=\mathbf{P}(V_{0},s_{0})$, is defined from which the spatio-
temporal perturbations are sampled for each $k$ as
$e_{k}=\mathbf{P}_{0}^{1/2}\zeta_{k}$, where $\zeta_{k}$ is a sample of a
centered and normalized Gaussian random vector, and where
$\mathbf{P}_{0}^{1/2}$ stands for the square root matrix of $\mathbf{P}_{0}$,
i.e.
$\mathbf{P}_{0}=\mathbf{P}_{0}^{1/2}\left(\mathbf{P}_{0}^{1/2}\right)^{\mathrm{T}}$.
The square root $\mathbf{P}_{0}^{1/2}$ has been computed from the singular
value decomposition of the matrix $\mathbf{P}_{0}$.

An example of a perturbation sample is presented in Fig. 2 where the temporal
evolution $e(t,x=0)$ is shown in panel (a) while the initial condition within
the domain, $e(t=0,x)$ is given in panel (b). Note that the time axis in panel
(a) has been inverted so to facilitate the understanding. The blue dots
corresponds to the value of the sampled error field $e$ at $t=0$ and $x=0$.

The figure Fig. 3 shows both variance and length-scale fields that are
computed from the PKF and the ensemble simulations and compared at three
different timestamps. The first panels (a) and (b) respectively show the
variance and the length-scale normalized by $L_{h}$ at initial time. As
prescribed, the initial variance is homogeneous and equal to $1$. The initial
length-scale is also homogeneous and equal to $L_{h}$. The panels (c) and (d)
present the evolution at time $t=0.25T_{adv}$. As expected the variance (the
length-scale) increases (decreases) according to the specified boundary
condition close to $x=0$ (red dot). The length-scale is also modified over the
whole domain, this is the effect of the velocity gradient of the term
$2s_{c,xx}\partial_{x}u$ in Eq. (21c). Since the velocity field is positive,
the variance and the length-scale are transported to the right of the domain.
Finally, panels (e) and (f) show the fields at time $t=1.5t_{adv}$. The
information injected by the boundary condition at $x=0$ has reached the other
side of the domain unscathed for the variance and length-scale fields.

In order to strengthen these results, we show in Fig. 4 the evolution through
time of the fields for the middle point in the domain $x=0.5\Lambda$. As
expected, the variance in panel (a) remains constant until the information
from the boundary condition arrives, where oscillations start, following the
prescribed sine shape of the boundary condition shifted in time. In panel (b),
the length-scale follows the same kind of dynamic except that the length-scale
varies from $t=0$ to $t=0.5T_{adv}$, a variation that is not due to the
boundary condition but to the heterogeneity of the wind field. Note that
ensemble estimation of the variance and of the length-scale are subject to
some sampling noise even with the large ensemble size $N_{e}=6400$.

Overall, this simulation shows no numerical artifact and the PKF and EnKF
forecasts overlap perfectly. Moreover, the continuous and differentiable error
statistics of the EnKF statistics shows that the generated duets of errors for
the initial condition and boundary condition have been appropriately
specified.

These results validate that the specification of the PKF boundaries proposed
in Section 3.1 is correct when applying Dirichlet condition in a transport
dynamics. Moreover it also validates the specification of the perturbations
Eq. (32), introduced in Section 4, for the ensemble validation to build
prescribed error statistics.

Note that, this example has also shown the ability of the PKF to apply for
open boundary condition.

Now, we validate the PKF boundary conditions applied for a diffusion equation.

###  5.3 Application to the diffusion equation

In this experiment setting, the heterogeneous diffusion Eq. (22) is
considered. The temporal discretization scheme used for the ensemble
simulation is a backward Euler scheme (implicit Euler method) with a fixed
time-step $dt_{BE}\approx 2.10^{-4}$. For the PKF dynamics we used a Runge-
Kutta scheme of order 4 with a fixed time-step $dt_{RK4}\approx 5.10^{-6}$.
The simulation is performed from time $t=0$ to $t_{end}=1.2T_{\mathrm{diff}}$
with $T_{\mathrm{diff}}=\frac{\Lambda^{2}}{4D_{max}}$ the time scale of the
diffusion of a half-domain.

The diffusion coefficient for the simulation is set as the heterogeneous
stationary field $D(x)=1+\frac{A}{A_{max}}$ with $A(x)=sin(\pi x)(1+x)^{8}$
where $A_{\mathrm{max}}=Max_{x}A(x)$, and is shown in Fig. 5. This diffusion
field reproduces the kind of diffusion encountered in the dynamics of
radiation belts in order to evaluate the ability of PKF to solve this problem.

Figure 5: Heterogeneous diffusion coefficient generated for the experiment

####  5.3.1 Dirichlet boundary conditions

To generate a coherent set of perturbations for the ensemble simulation, the
same technique described in Section 5.2 is considered except that both
boundaries at $x=0$ and $x=\Lambda$ are subject to Dirichlet conditions. The
extended domain considered is
$[-\sqrt{D(0)t_{end}},0]\cup[0,\Lambda]\cup[\Lambda,\Lambda+\sqrt{D(\Lambda)t_{end}}]$.

This time, the parameters considered for the simulation and ensemble
generation are as follows, the initial variance is set to the linear function
$V_{0}(t=0,x)=1+\frac{3}{\Lambda}x$ and the initial length-scale is set
homogeneous and equal to 10% of the domain length
$L_{0}(t=0,x)=L_{h}=0.1\Lambda$. For the left boundary condition at $x=0$, the
variance and the length-scale are stationary and set equal to $1$ and $L_{h}$
respectively i.e. $V_{0}(t,x=0)=1$ and $L_{0}(t,x=\Lambda)=L_{h}$. For the
right boundary condition at $x=\Lambda$, the variance and the length-scale are
stationary and set equal to $4$ and $L_{h}$ respectively i.e.
$V_{0}(t,x=\Lambda)=4$ and $L_{0}(t,x=\Lambda)=L_{h}$. From this
specification, an ensemble of perturbations has been populated following the
same procedure, $e_{k}=\mathbf{P}_{0}^{1/2}\zeta_{k}$, as detailed in Section
5.2. The resulting perturbations are similar to the ones shown in Fig. 2 for
the advection, except that there is a right extension of the domain in
addition of the left extension for the advection (not shown).

Figure 6: Comparison of the forecast-error variance (left column) and
normalized length-scale (right column) fields dynamics for the heterogeneous
diffusion equation on a 1D bounded domain with Dirichlet boundary conditions,
and shown at times $t=0$, $t=0.2T_{\mathrm{diff}}$ and
$t=1.5T_{\mathrm{diff}}$.  Figure 7: Time evolution of the forecast-error
variance (a) and normalized length-scale (b) at $x=0.5\Lambda$, for the
diffusion equation with Dirichlet boundary conditions.

The comparison between the PKF and EnKF predictions at different time steps
are shown in Fig. 6. The first panels (a) and (b) are coherent with the
specification of the initial condition for both the EnKF and the PKF. Panels
(c) and (d) show the evolution of the variance and length-scale at
$t=0.2T_{\mathrm{diff}}$.

Due to physical diffusion, far from the boundaries e.g. at the center of the
domain, the magnitude of the error is expected to decrease over time with an
attenuation of the variance, while the length-scale should increase; and at
the boundaries the uncertainty should remained as specified by the Dirichlet
conditions. This is precisely the behaviour observed for both the EnKF and the
PKF, at the center of the domain and at the boundaries where the Dirichlet
condition imposes fixed values for the variance and the length-scale on both
sides of the domain.

However, panel (d) shows a noticeable gap between the length-scale computed by
the PKF and the one estimated from the ensemble. This gap can be due to the
closure Eq. (23) but it has a limited impact on the variance field (panel c)
which suggests that the PKF prediction of the variance is an accurate proxy
for the EnKF estimation.

On the last panels (e) and (f), the variance and the length-scale settle down
and the values predicted by the PKF are close to the values computed from the
ensemble except for the error observed between the length-scale fields in the
middle of the domain. As seen in Fig. 7, the variance and length-scale are
close to the permanent regime at $t=1.2T_{\mathrm{diff}}$ showing that the PKF
performed well even over a significant time period.

To conclude, this experiment has confirmed the specification of the Dirichlet
boundary conditions of Section 3.1 for the PKF applied to a heterogeneous
diffusion equation. It has shown the ability of the PKF to accurately
approximate the uncertainty dynamics as diagnosed from the EnKF but at a lower
cost corresponding the price of two time integrations compared to the $6400$
integrations needed for the ensemble. Another result is that the simulations
also validate the theoretical derivation of the time-scale setting Eq. (33)
needed to obtain a specific length-scale at the boundaries.

We end the numerical validation by considering the Neumann conditions applied
to the heterogeneous diffusion equation.

####  5.3.2 Neumann boundary conditions

Figure 8: Samples of random error generated as initial condition that verify
the Neumann condition at the boundaries $x=0$ and $x=\Lambda$.  Figure 9:
Forecast-error metric field for the heterogeneous diffusion equation on a 1D
bounded domain with Neumann boundary conditions, shown at times $t=0$,
$t=0.2T_{\mathrm{diff}}$ and $t=1.2T_{\mathrm{diff}}$.  Figure 10: Comparison
of the forecast-error variance (left column) and length-scale (right column)
fields dynamics for the heterogeneous diffusion equation on a 1D bounded
domain with Neumann boundary conditions, and shown at times $t=0$,
$t=0.2T_{\mathrm{diff}}$ and $t=1.2T_{\mathrm{diff}}$.  Figure 11: Time
evolution of the forecast-error variance (a) and length-scale (b) at
$x=0.5\Lambda$, for the diffusion equation with Neumann boundary conditions.

As above mentioned in Section 4, compared with the Dirichlet, the Neumann
conditions are simulated in an ensemble of forecasts, as an initial condition
problem without perturbation at the boundaries. The problem is then to produce
an ensemble of initial conditions that verify the Neumann conditions.

To do so, a covariance model based on a homogeneous pseudo-diffusion equation
has been considered [52]. The terminology pseudo means that the diffusion is
not physical but only a tricky way to create large covariance model as used in
variational data assimilation. In particular the square-root covariance
$\mathbf{P}_{0}^{1/2}$ resulting from the integration of the pseudo-diffusion
equation reads as the linear operator

$\mathbf{P}_{0}^{1/2}=\mathbf{\Sigma WL},$ (34)

where $\mathbf{L}=e^{\frac{1}{2}\kappa\partial_{x}^{2}}$ is the propagator
associated with the diffusion equation

$\partial_{\tau}u=\kappa\partial_{x}^{2}u$ (35)

of pseudo-time $\tau$, integrated from $\tau=0$ to $\tau=\frac{1}{2}$, and
using Neumann conditions at the boundaries [26]; $\mathbf{W}$ is a diagonal
normalisation so that $\mathbf{WL}(\mathbf{WL})^{\mathrm{T}}$ is a correlation
operator ; and $\mathbf{\Sigma}$ is a diagonal matrix of standard deviations,
so that the spatial variance field is the linear profile with
$V_{0}(t=0,x=0)=1$ and $V_{0}(t=0,x=\Lambda)=4$. Note that the pseudo-
diffusion coefficient $\kappa$ is related to the length-scale $l$ of the
correlation functions as to $\kappa=l^{2}/2$ [33]. For the numerical
application, $l_{h}=0.1\Lambda$.

Again, an ensemble of initial conditions are populated from the square-root
Eq. (34), $e_{k}=\mathbf{P}_{0}^{1/2}\zeta_{k}$, where $\zeta_{k}$ is a sample
of centered Gaussian random vector. Fig. 8 shows some samples of the
normalized error resulting from Eq. (34) i.e.
$\varepsilon_{k}=\mathbf{WL}\zeta_{k}$. As it is expected, the normalized
error are flat at the boundaries (red arrows pointing toward the interior of
the domain). The resulting anisotropy diagnosed from the ensemble of initial
condition $t=0$ leads to the metric field shown in Fig. 9 (in blue but super-
imposed by the orange line). As expected, far from the boundary, the metric is
homogeneous equal to $g=1/l_{h}^{2}$ i.e. near $x=0.5\Lambda$, but oscillates
near the boundaries to reach a value of zero at the boundaries. The
oscillations is due to the constraint of symmetry of the covariance matrix
[32].

As discussed in Section 3.2, for Neumann conditions the PKF dynamics is solved
following its metric formulation, which is given by Eq. (25) for the physical
diffusion equation Eq. (22). For the numerical validation of the Neumann BCs,
the initial condition for the PKF is the variance field of linear profile
shown in Fig. 10-(a) and the experiment metric field diagnosed from the
ensemble of initial conditions shown in Fig. 9 for $t=0$ (orange line,
superimposed to the blue line of the EnKF diagnosis).

The PKF dynamics is computed and the results are compared with the ensemble of
forecasts of the heterogeneous diffusion equation Eq. (22) and Neumann
conditions on both sides. The results are shown in Fig. 9 for the variance and
the length-scale (computed from the inverse of the metric), and in Fig. 9 for
the metric. The results are shown for times of interest selected from the time
evolution reproduced in Fig. 11 where a relaxation toward a stationary state
of uncertainty appears.

As expected for a diffusion, the variance decreases along the time, while the
length-scale increases. Note that for Neumann condition, the variance at the
boundary also decreases while it was constant in the Dirichlet condition. For
the ensemble estimation, the length-scales at the boundaries (blue lines in
panel (b-d-f) ) are large but finite where it is expected to be infinite: this
is due to the numerical estimation of the length-scale deduced from Eq. (14),
while the metric remains null at the boundaries during the simulation (see
Fig. 9 the red dots).

Compared with the EnKF diagnosis, the PKF performs well by reproducing the
same behaviour of the uncertainty dynamics as for the EnKF, except that the
length-scale predicted from the PKF underestimates the length-scale diagnosed
from the ensemble. However, the very large length-scale values, larger than
the domain size $\Lambda$, as diagnosed from the ensemble is subject to the
limitation of the numerical computation of Eq. (14) for large correlations
that can present a positive bias [31]. Moreover, the large length-scale of the
EnKF can also be influenced by model error [34]. Because of these limitations,
it is not certain that the EnKF reproduces the true dynamics of the
uncertainty for these extreme values of the length-scales, while it is
considered as the reference. Hence, the discrepancy between the PKF and the
EnKF reference, may not be due to a defect of the PKF that could be better
than the ensemble estimation here.

To conclude, this experiment has confirmed the specification of the Neumann
boundary conditions of Section 3.2 for the PKF applied to a heterogeneous
diffusion equation. It has shown the ability of the PKF to accurately
approximate the uncertainty dynamics as diagnosed from the EnKF but at a lower
cost corresponding the price of two time integrations compared to the $6400$
integrations needed for the ensemble.

This ends the validation of the specification of the boundary conditions for
the PKF. The summary of the results obtained in the paper as well as the
perspective of the work are given in the following section of conclusion.

##  6 Conclusion

This work contributed to explore the parametric Kalman filter (PKF), that is a
recent approximation of the Kalman filter proposed for application in large
systems. The parametric approach investigated here consists to approximate the
forecast-error covariance matrix by a covariance model parameterized from the
variance and the anisotropy. The anisotropy can be specified in term of metric
tensor or its inverse, the aspect tensor, that is the square of the length-
scale in 1D domains. The PKF dynamics describes how the mean, the variance and
the anisotropy evolve in time, leading to the prediction of the error
statistics at a lower numerical cost than the one needed for the full
covariance propagation described in the Kalman filter forecast step.

In this contribution, we proposed how to specify the error statistics at the
boundary of a domain when considering a PKF forecast of the uncertainty. We
detail here pragmatic solutions for large systems with strong variability at
their domain’s edge, such as in air quality at regional scale or in radiation
belts "weather".

Two kind of boundaries have been considered, the Dirichlet and the Neumann
conditions depending on the dynamics. We obtained that a Dirichlet boundary
condition for the dynamics implies to specify a Dirichlet conditions for the
variance and the anisotropy in the PKF dynamics, no mater how the PKF dynamics
is specified in metric or in aspect tensor. Similarly, a Neumann condition
implies to specify a Neumann condition for the variance but a Dirichlet
condition for the metric. However, for Neuman conditions, the formulation of
the PKF in metric is more adapted than the one in aspect tensor, which would
require infinite value at the boundary.

Two dynamics of interest for weather forecasting, air quality or radiation
belt dynamics have been considered: the transport and the diffusion equation.
For these dynamics, the PKF predictions have been validated numerically in a
1D limited area domain, from the comparison to an ensemble estimation e.g. by
considering heterogeneous and non-stationnary uncertainty at the boundary for
Dirichlet conditions. For Dirichlet conditions, the ensemble of forecast has
been designed by adding, to each member, a perturbation of the boundary with a
varying correlation time-scale, for which we derived analytical formula that
depend on the dynamics and introduced to obtain a prescribed spatial
anisotropy near the boundary. For both the advection and the diffusion, the
PKF has been shown able to reproduce the uncertainty dynamics diagnosed from
the ensemble of forecast. Note that this result is not obvious, even when
considering a linear-Gaussian setting: the PKF based on VLATcov model only
approximate the full covariance dynamics from the dynamics of the variance and
the anisotropy, that is a crude approximation of the full covariance dynamics.
Hence, it is an encouraging result for the application in air quality (for the
transport equation) and in radiation belts monitoring (for the diffusion). For
other applications, these results will have to be confirmed by considering
non-linear dynamics of interest, and multivariate statistics extension of the
PKF, as it is required for numerical weather prediction ; which needs
complementary exploration of the PKF beyond the scope of the present
contribution and are challenging topics. Note also that the numerical testbed
indirectly validates the auto-correlation time scale formula, and constitutes
a contribution to the ensemble methods (while these formula are is not needed
for the PKF).

The next step will be to study the BC conditions for domains of larger
dimensions, where we expect some changes e.g. non-zero components of the
metric tensor along the tangential direction to the boundary in Neumann
conditions.

We can mention that the dynamics of the uncertainty for bounded domains can be
of importance in variational data assimilation or observation targeting
applied for local area models, that could be another topic to investigate with
the PKF.

Beyond these challenging topics, we can mention that the results in 1D should
already find important applications e.g. in the dynamics of uncertainty in the
boundary layer for air quality, or wild-land fire predictions ; or in the
exploration of the coupling of uncertainty in the coupled atmosphere-ocean
system during the alternating integration of the two fluids.

##  Appendix A Accounting for the uncertainty at the boundary condition in a
KF

This section illustrates the KF equation in presence of a boundary, applied to
a transport. The KF forecast step is first introduced where the boundary
condition appears as a command. Then it is applied when the boundary is
perfectly known or uncertain. An illustration of uncertain boundary is
proposed when the error is a stationary auto-regressive noise of order 1 (
AR(1) ), for which the formalism can be reworded without command but
considering an extended domain with an equivalent model error.

###  A.1 KF forecast step for a transport in a limited-area 1D domain

The transport $\partial_{t}c+U\partial_{x}c=0$, of a concentration field
$c(t,x)$ at a constant positive velocity $U$, along a 1D limited-area domain
$[0,\Lambda]$ of coordinate $x$, with a Dirichlet condition $c(t,x=0)=b(t)$ at
the inlet and open condition at the outlet $x=1$, discretized on a regular
grid in time $(t_{q})_{q\in\mathbb{N}}$ and space $(x_{k})_{k\in[0,n-1]}$
(i.e. $x_{0}=0$) with a time step $\delta t$ and a grid space $\delta x$
defined by $\delta t=\frac{\delta x}{U}$, and when the boundary condition is
constant over a time step, the dynamics reads as

$c_{q+1}=\mathbf{M}c_{q}+\mathbf{F}b_{q},$ (36)

where here, $c_{q+1}$ denote the vector resulting from the discretization of
$c(t_{q},x)$, $\mathbf{F}=[1,0,\cdots,0]^{T}$ and

$\mathbf{M}=\begin{bmatrix}0&0&0&\cdots&0\\\ 1&0&0&\cdots&0\\\
0&1&0&\cdots&0\\\ \vdots&\ddots&\ddots&\ddots&\vdots\\\
&\cdots&0&1&0\end{bmatrix},$ (37)

so that the value of the concentration $c(t_{q+1},x_{k})$ at a grid point
$x_{k}$ and at a time $t_{q+1}$ is equals to $c(t_{q},x_{k-1})$ providing that
$0<k<n-1$ ; at the boundaries the the first component of $c_{q+1}$ is $b_{q}$
(i.e. $c(t_{q+1},0)=b_{q}$) while the value of $c(t_{q},x_{n-1})$. After a
long period of time integration, when $t_{q}>\Lambda/U$,
$c(t_{q+1},x_{k})=b_{q-k}$ for any $k\in[0,n-1]$. Note that, in optimal
control, $b_{q}$ in Eq. (36) plays the role of a command.

From Eq. (36), it follows that the dynamics of the forecast error is written
as

${e^{f}_{c}}_{q+1}=\mathbf{M}{e^{f}_{c}}_{q}+\mathbf{F}{e_{b}}_{q},$ (38)

where $e^{f}_{c}$ ($e_{b}$) denotes the forecast error on $c$ (the error on
$b$). As it is common in KF theory, the errors are assumed centered i.e.
$\mathbb{E}\left[e^{f}_{c}\right]=\mathbb{E}\left[e_{b}\right]=0$ for any
time).

###  A.2 Perfectly known boundary condition

In the case where the boundary condition at the inlet is perfectly known (i.e.
${e_{b}}_{q}=0$ for any $q$), but the initial condition is uncertain,
characterized by Gaussian error distribution of covariance $\mathbf{P}^{a}$ ;
then the Gaussian uncertainty, characterized by the covariance
$\mathbf{P}^{f}_{q}$ at $t_{q}$, evolves following the Kalman filter forecast
step Eq. (4) that is

$\mathbf{P}^{f}_{q+1}=\mathbf{M}\mathbf{P}^{f}_{q}\mathbf{M}^{T}.$ (39)

Hence, after one time step ($q=1$), the forecast-error covariance matrix reads
as $\mathbf{P}^{f}_{1}=\mathbf{M}\mathbf{P}^{a}\mathbf{M}^{T}$
($\mathbf{P}^{f}_{0}=\mathbf{P}^{a}$), that is a matrix that looks like

$\mathbf{P}^{f}_{1}=\begin{bmatrix}0&0&0&\cdots&0\\\ 0&&&&\\\
\vdots&&{(\mathbf{P}^{f}_{0})}^{-}&&\\\ 0&&&&\end{bmatrix}$ (40)

where ${(\mathbf{P}^{f}_{0})}^{-}$ denotes square matrix of dimension $n-1$
whose coefficients are those of the submatrix $\mathbf{P}^{a}_{[0:n-2,0:n-2]}$
of $\mathbf{P}^{a}$ since $\mathbf{P}^{f}_{0}=\mathbf{P}^{a}$. It appears that
the variance at the boundary is null, which is coherent when the boundary
condition (not necessary constant in time) is perfectly known. After, $q>0$
iterations of Eq. (39),
$\mathbf{P}^{f}_{q+1}=(\mathbf{M})^{q}\mathbf{P}^{f}_{0}(\mathbf{M}^{T})^{q}$
, the $q^{th}$ first rows and lines of $\mathbf{P}^{f}_{q}$ are null. Since
$\mathbf{M}$ is nilpotent of order $n$ (i.e. $\mathbf{M}^{n-1}\neq 0$ and
$\mathbf{M}^{n}=0$), after $n$ iterations, the forecast-error covariance
matrix is null, meaning that there is no more error within the domain
(independently of the initial uncertainty of the initial condition).

While in this section we considered a perfectly known boundary, this is not
the situation encountered in real applications e.g. in forecasting at regional
scale where uncertainty at large scale should introduce an inflow of
uncertainty at the entrance of the domain.

Hence, the inflow of "certainty" that results from the common covariance
propagation Eq. (39) of the Kalman filter forecast step, but used without
care, explains the loss of variance observed in ensemble forecasting when no
perturbation of the boundary is used to represent the uncertainty [27, 28].

The next section provides the correct formulation of the covariance dynamics
in presence of uncertainty at the boundary.

###  A.3 Uncertain boundary condition

Now, when considering that the boundary condition is uncertain, the forecast-
error covariance matrix evolves as

$\mathbf{P}^{f}_{q}=\mathbf{M}\mathbf{P}^{f}_{q}\mathbf{M}^{T}+\mathbf{F}\mathbb{E}\left[{e_{b}}_{q}({e_{b}}_{q})^{T}\right]\mathbf{F}^{T}+\mathbf{M}\mathbb{E}\left[{e_{c}^{f}}_{q}({e_{b}}_{q})^{T}\right]\mathbf{F}^{T}+\left(\mathbf{M}\mathbb{E}\left[{e_{c}^{f}}_{q}({e_{b}}_{q})^{T}\right]\mathbf{F}^{T}\right)^{T},$
(41)

and reads as

$\mathbf{P}^{f}_{q}=\begin{bmatrix}0&0&\cdots&0\\\ 0&&&\\\
\vdots&&{(\mathbf{P}^{f}_{q})}^{-}&\\\
0&&&\end{bmatrix}+\begin{bmatrix}{V_{b}}_{q}&0&\cdots&0\\\ 0&&&\\\
\vdots&&{(0)}&\\\ 0&&&\end{bmatrix}+\begin{bmatrix}0&0&\cdots&0\\\ *&&&\\\
\vdots&&{(0)}&\\\ *&&&\end{bmatrix}+\begin{bmatrix}0&*&\cdots&*\\\ 0&&&\\\
\vdots&&{(0)}&\\\ 0&&&\end{bmatrix},$ (42)

where ${V_{b}}_{q}=\mathbb{E}\left[{e_{b}}_{q}({e_{b}}_{q})^{T}\right]$ is the
variance of the uncertainty at the boundary, and where stars, $*$, correspond
to a priori non zero coefficients which represents the correlations between
the uncertainty at the boundary with the forecast-error in the remaining part
of the domain $(0,\Lambda]$. Hence, the first term in Eq. (45),
$\mathbf{M}\mathbf{P}^{f}_{q}\mathbf{M}^{T}$, is similar to the matrix shown
in Eq. (40). The other terms in Eq. (45) will contribute to fill the first row
and column of the matrix: the second term of Eq. (45),
$\mathbf{F}\mathbb{E}\left[{e_{b}}_{q}({e_{b}}_{q})^{T}\right]\mathbf{F}^{T}={V_{b}}_{q}\mathbf{F}\mathbf{F}^{T}$,
where $\mathbf{F}\mathbf{F}^{T}$ is the null matrix except for the first
coefficient that is one, fills the first coefficient of
$\mathbf{P}^{f}_{q+1}$. Similarly, the third and the fourth terms of Eq. (45)
will fill the first row and the first column.

###  A.4 Example with a stationary AR(1) boundary error

As an instructive example, we consider the case where the error is a
stationary AR(1) random process,
${e_{b}}_{q+1}=a{e_{b}}_{q}+\sqrt{1-a^{2}}\sqrt{V_{b}}\zeta_{q}$, where
$V_{b}$ denotes the time independent variance of ${e_{b}}_{q}$ (independent of
$q$) ; $a=\exp{\left(-\delta t/\tau\right)}$ where $\tau$ is the correlation
time scale ; and $\zeta=(\zeta_{q})_{q\in\mathcal{N}}$ is a process of
independent a nd identical independent random variable, where for each $q$,
$\zeta_{q}$ is normal law i.e. $\zeta_{q}\sim\mathcal{N}(0,1)$. The time
correlation for this processes is the exponential
$cor({e_{b}}_{q},{e_{b}}_{q+r})=\exp{\left(-r\delta t/\tau\right)}$, that also
reads as $cor({e_{b}}_{q},{e_{b}}_{q-r})=\exp{\left(-r\delta t/\tau\right)}$
because of the stationarity of the process (homogeneity in time).

Since after a long period, $c(t_{q+1},x_{k})=b_{q-k}$, it follows that
${e_{c}}(t_{q+1},x_{k})={e_{b}}_{q-k}$ so that the $k^{th}$ component of the
vector $\mathbb{E}\left[{e_{c}^{f}}_{q}({e_{b}}_{q})^{T}\right]$ is
$cov(e_{c}^{f}(t_{q},x_{k}),{e_{b}}_{q})=V_{b}cor({e_{b}}_{q},{e_{b}}_{q-k})=V_{b}cor({e_{b}}_{q+k},{e_{b}}_{q})=V_{b}cor({e_{b}}_{q+k},{e_{b}}_{q})=V_{b}\exp{\left(-k\delta
t/\tau\right)}$. With $e_{c}^{f}(t_{q+1},0)={e_{b}}_{q}$,
$e_{c}^{f}(t_{q+1},x_{k})=e_{c}^{f}(t_{q},x_{k-1})$, and $\delta t=\delta
x/U$, then
$cov(e_{c}^{f}(t_{q+1},x_{0}),e_{c}^{f}(t_{q+1},x_{k}))=V_{b}\exp{\left(-k\delta
x/(U\tau)\right)}$. Hence, after a relaxation period, the error covariance
matrix is stationnary with a spatial correlation function given as an
exponential correlation of length-scale $l=U\tau$.

###  A.5 Extension of the space domain for a boundary error modeled by an
AR(1)

When the boundary error is an auto-regressive process or order $1$, it is
possible to rewords the dynamics Eq. (38), by extending the physical domain
with an additional one point corresponding to the boundary condition that will
enter to the domain in the future. Focusing on error dynamics, this reads as
follow:

Considering the extended vector $\tilde{e}_{q}=[{e_{b}}_{q},e^{f}_{q}]^{T}$ of
dimension $n+1$, then Eq. (38) is written as

$\tilde{e}_{q+1}=\widetilde{\mathbf{M}}\tilde{e}_{q}+\tilde{\mathbf{F}}\zeta_{q},$
(43)

where

$\widetilde{\mathbf{M}}=\begin{bmatrix}a&0\\\
\mathbf{F}&\mathbf{M}\end{bmatrix},$ (44)

and $\widetilde{\mathbf{F}}=[\sqrt{1-a^{2}}V_{b},0,\cdots,0]^{T}$. With this
formulation, Eq. (43) is similar to the classic formulation of the KF equation
in presence of model error (that would be here the second term
$\tilde{\mathbf{F}}\zeta_{q}$), leading to the error covariance dynamics

$\widetilde{\mathbf{P}}^{f}_{q}=\widetilde{\mathbf{M}}\widetilde{\mathbf{P}}^{f}_{q}\widetilde{\mathbf{M}}^{T}+\widetilde{\mathbf{Q}},$
(45)

where
$\widetilde{\mathbf{Q}}=\widetilde{\mathbf{F}}\widetilde{\mathbf{F}}^{T}$.

This extension of the domain also applies for auto-regressive process of order
$m$ where we need to add $m$ points.

##  Appendix B Closure of the PKF Dynamics for the diffusion equation

The computation of the PKF dynamics for the diffusion equation Eq. (22), with
SymPKF, leads to the dynamical system

$\displaystyle\partial_{t}f$
$\displaystyle=D\partial^{2}_{x}f+\partial_{x}D\partial_{x}f,$ (46a)
$\displaystyle\partial_{t}V_{f}$
$\displaystyle=-\frac{2DV_{f}}{s_{f,xx}}+D\partial^{2}_{x}V_{f}-\frac{D\left(\partial_{x}V_{f}\right)^{2}}{2V_{f}}+\partial_{x}D\partial_{x}V_{f},$
(46b)
$\partial_{t}s_{f,xx}=2Ds_{f,xx}^{2}{\mathbb{E}}\left(\varepsilon_{f}\partial^{4}_{x}\varepsilon_{f}\right)-3D\partial^{2}_{x}s_{f,xx}-2D+\\\
\frac{6D\left(\partial_{x}s_{f,xx}\right)^{2}}{s_{f,xx}}-\frac{2Ds_{f,xx}\partial^{2}_{x}V_{f}}{V_{f}}+\frac{D\partial_{x}V_{f}\partial_{x}s_{f,xx}}{V_{f}}+\\\
\frac{2Ds_{f,xx}\left(\partial_{x}V_{f}\right)^{2}}{V_{f}^{2}}-2s_{f,xx}\frac{d^{2}}{dx^{2}}D+\\\
2\partial_{x}D\partial_{x}s_{f,xx}-\frac{2s_{f,xx}\partial_{x}D\partial_{x}V_{f}}{V_{f}}$
(46c)

where this time the term
$\mathbb{E}\left[\varepsilon_{f}\partial^{4}_{x}\varepsilon_{f}\right]$ is not
determined from $f$, $V_{f}$ and $s_{f,xx}$. This dynamics can be closed
considering the closure Eq. (23).

##  Appendix C Specification of the temporal metric tensor for evolution
equations

This section details the link between the temporal metric Eq. (30a),
$\mathbf{g}_{tt}=\mathbb{E}\left[\partial_{t}\varepsilon\partial_{t}\varepsilon\right]$,
and the dynamics of the error. Since the trend of the normalized error reads
as

$\partial_{t}\varepsilon=\frac{1}{\sqrt{V}}\partial_{t}e-\frac{1}{2V^{3/2}}e\partial_{t}V,$
(47)

then the temporal metric tensor is written as

$g_{tt}=\frac{1}{V}\mathbb{E}\left[\left(\partial_{t}e\right)^{2}\right]-\frac{1}{V^{2}}\mathbb{E}\left[e\partial_{t}e\right]\partial_{t}V+\frac{1}{4V^{3}}\mathbb{E}\left[e^{2}\right]\left(\partial_{t}V\right)^{2}.$
(48)

However, we recognize the expression of the variance
$V=\mathbb{E}\left[e^{2}\right]$ and its trend, Eq. (16a), so that the
temporal metric simplifies as

$g_{tt}=\frac{1}{V}\mathbb{E}\left[\left(\partial_{t}e\right)^{2}\right]-\frac{1}{4V^{2}}\left(\partial_{t}V\right)^{2}.$
(49a) Introducing the trend of the error Eq. (3) and by definition of
$\varepsilon=e/\sqrt{V}$, the temporal metric reads as
$g_{tt}=\frac{1}{V}\mathbb{E}\left[\left(\mathcal{M}(\varepsilon\sqrt{V},\partial(\varepsilon\sqrt{V}))\right)^{2}\right]-\frac{1}{4V^{2}}\left(\partial_{t}V\right)^{2}.$
(49b)

##  Appendix D Time auto-correlation boundary condition for the diffusion
equation

The computation of the time auto-correlation metric Leverages on SymPKF. For
the diffusion equation, SymPKF leads to

$g_{f,tt}=D^{2}{\mathbb{E}}\left(\varepsilon_{f}\partial^{4}_{x}\varepsilon_{f}\right)+2D^{2}\partial^{2}_{x}g_{f,xx}-\frac{D^{2}g_{f,xx}\partial^{2}_{x}V_{f}}{V_{f}}+\\\
\frac{D^{2}\partial_{x}V_{f}\partial_{x}g_{f,xx}}{V_{f}}+\frac{3D^{2}g_{f,xx}\left(\partial_{x}V_{f}\right)^{2}}{2V_{f}^{2}}+\frac{D^{2}\left(\partial^{2}_{x}V_{f}\right)^{2}}{4V_{f}^{2}}\\\
-\frac{D^{2}\left(\partial_{x}V_{f}\right)^{2}\partial^{2}_{x}V_{f}}{4V_{f}^{3}}+\frac{D^{2}\left(\partial_{x}V_{f}\right)^{4}}{16V_{f}^{4}}+D\partial_{x}D\partial_{x}g_{f,xx}+\\\
\frac{Dg_{f,xx}\partial_{x}D\partial_{x}V_{f}}{V_{f}}+\frac{Dg_{f,xx}\partial_{t}V_{f}}{V_{f}}+\frac{D\partial_{x}D\partial_{x}V_{f}\partial^{2}_{x}V_{f}}{2V_{f}^{2}}\\\
-\frac{D\partial_{t}V_{f}\partial^{2}_{x}V_{f}}{2V_{f}^{2}}-\frac{D\partial_{x}D\left(\partial_{x}V_{f}\right)^{3}}{4V_{f}^{3}}+\frac{D\partial_{t}V_{f}\left(\partial_{x}V_{f}\right)^{2}}{4V_{f}^{3}}+\\\
g_{f,xx}\left(\partial_{x}D\right)^{2}+\frac{\left(\partial_{x}D\right)^{2}\left(\partial_{x}V_{f}\right)^{2}}{4V_{f}^{2}}-\frac{\partial_{x}D\partial_{t}V_{f}\partial_{x}V_{f}}{2V_{f}^{2}}+\\\
\frac{\left(\partial_{t}V_{f}\right)^{2}}{4V_{f}^{2}}.$ (50)

Considering the analytical closure Eq. (23) for the unclosed term
${\mathbb{E}}\left(\varepsilon_{f}\partial^{4}_{x}\varepsilon_{f}\right)$, the
correspondence is written as

$g_{f,tt}=3D^{2}g_{f,xx}^{2}-\frac{D^{2}g_{f,xx}\partial^{2}_{x}V_{f}}{V_{f}}+\frac{D^{2}\partial_{x}V_{f}\partial_{x}g_{f,xx}}{V_{f}}+\\\
\frac{3D^{2}g_{f,xx}\left(\partial_{x}V_{f}\right)^{2}}{2V_{f}^{2}}+\frac{D^{2}\left(\partial^{2}_{x}V_{f}\right)^{2}}{4V_{f}^{2}}-\frac{D^{2}\left(\partial_{x}V_{f}\right)^{2}\partial^{2}_{x}V_{f}}{4V_{f}^{3}}+\\\
\frac{D^{2}\left(\partial_{x}V_{f}\right)^{4}}{16V_{f}^{4}}+D\partial_{x}D\partial_{x}g_{f,xx}+\frac{Dg_{f,xx}\partial_{x}D\partial_{x}V_{f}}{V_{f}}+\\\
\frac{Dg_{f,xx}\partial_{t}V_{f}}{V_{f}}+\frac{D\partial_{x}D\partial_{x}V_{f}\partial^{2}_{x}V_{f}}{2V_{f}^{2}}\\\
-\frac{D\partial_{t}V_{f}\partial^{2}_{x}V_{f}}{2V_{f}^{2}}-\frac{D\partial_{x}D\left(\partial_{x}V_{f}\right)^{3}}{4V_{f}^{3}}+\frac{D\partial_{t}V_{f}\left(\partial_{x}V_{f}\right)^{2}}{4V_{f}^{3}}+\\\
g_{f,xx}\left(\partial_{x}D\right)^{2}+\frac{\left(\partial_{x}D\right)^{2}\left(\partial_{x}V_{f}\right)^{2}}{4V_{f}^{2}}-\frac{\partial_{x}D\partial_{t}V_{f}\partial_{x}V_{f}}{2V_{f}^{2}}+\frac{\left(\partial_{t}V_{f}\right)^{2}}{4V_{f}^{2}}.$
(51)

The latter expression being quite complex, simplifications are introduced.
First the variance field is assumed locally homogeneous at the boundary i.e.
$\partial_{x}V_{f}(t,x=0)=0$, so that Eq. (51) simplifies as

$g_{f,tt}=3D^{2}g_{f,xx}^{2}+D\partial_{x}D\partial_{x}g_{f,xx}+\frac{Dg_{f,xx}\partial_{t}V_{f}}{V_{f}}+\\\
g_{f,xx}\left(\partial_{x}D\right)^{2}+\frac{\left(\partial_{t}V_{f}\right)^{2}}{4V_{f}^{2}}.$
(52)

Then, if the variance is moreover assumed stationary, then Eq. (52) becomes

$g_{f,tt}=3D^{2}g_{f,xx}^{2}+D\partial_{x}D\partial_{x}g_{f,xx}+g_{f,xx}\left(\partial_{x}D\right)^{2}.$
(53)

Eventually, then the diffusion coefficient field is homogeneous, then the
spatio-temporal connection between the temporal metric and the spatial metric
reads

$g_{f,tt}=3D^{2}g_{f,xx}.$ (54)

While Eq. (54) is a particular case, this equality is considered as a proxy
for setting the auto-correlation time scale of the boundary perturbation even
when the variance and the diffusion fields are heterogeneous.

Note that another expression for the spatio-temporal consistency Eq. (51) can
be obtained when first considering the dynamics of the variance given by Eq.
(24b), leading to replace the trend of the variance by
$\partial_{t}V_{f}=-2DV_{f}g_{f,xx}+D\partial^{2}_{x}V_{f}-\frac{D\left(\partial_{x}V_{f}\right)^{2}}{2V_{f}}+\partial_{x}D\partial_{x}V_{f}$,
so that Eq. (51) simplifies as

$g_{f,tt}=2D^{2}g_{f,xx}^{2}+\frac{D^{2}\partial_{x}V_{f}\partial_{x}g_{f,xx}}{V_{f}}+\\\
\frac{D^{2}g_{f,xx}\left(\partial_{x}V_{f}\right)^{2}}{V_{f}^{2}}+D\partial_{x}D\partial_{x}g_{f,xx}+\\\
\frac{2Dg_{f,xx}\partial_{x}D\partial_{x}V_{f}}{V_{f}}+g_{f,xx}\left(\partial_{x}D\right)^{2},$
(55)

from which the assumption of local homogeneity at the boundary i.e.
$\partial_{x}V_{f}(t,x=0)=0$, leads to

$g_{f,tt}=2D^{2}g_{f,xx}^{2}+D\partial_{x}D\partial_{x}g_{f,xx}+g_{f,xx}\left(\partial_{x}D\right)^{2}.$
(56)

When the diffusion field is constant, then the time auto-correlation metric is
related to the space auto-correlation metric by

$g_{f,tt}=2D^{2}g_{f,xx},$ (57)

which is a different result from Eq. (54).

It is not clear whether the appropriate consistency should be given by Eq.
(54) or Eq. (57) i.e. if it is right to replace the trend of the variance Eq.
(24b) in the consistency relation Eq. (51).

From numerical experiment, it appears that setting the time auto-correlation
of boundary perturbation with Eq. (54) in the EnKF is in agreement with the
PKF results. This suggests that taking into account the trend of the variance
would lead to a kind of over-specification of the boundary condition for the
diffusion equation in an EnKF approach.

##  Appendix E Open Research

V1.0 of the _Boundary conditions for the parametric kalman filter forecast_
software used to compute and analyze the numerical experiments presented in
this paper is preserved at 10.5281/zenodo.7193985 and developed openly at
https://github.com/opannekoucke/pkf-boundary. [41]

## Acknowledgements

The authors would like to thank François Rogier (ONERA) for fruitful
discussions. M. Sabathier PhD Thesis is supported by CNES grant 51/19168. O.
Pannekoucke is supported by the French national program LEFE/INSU grant
"Multivariate Parametric Kalman Filter" (MPKF).

## References

  * [1] Laaziz El Amraoui, Bojan Sič, Andrea Piacentini, Virginie Marécal, Nicolas Frebourg, and Jean-Luc Attié.  Aerosol data assimilation in the MOCAGE chemical transport model during the TRAQA/ChArMEx campaign: lidar observations.  Atmospheric Measurement Techniques, 13(9):4645–4667, sep 2020.
  * [2] C. Bishop, B. Etherton, and S. Majumdar.  Adaptive sampling with the ensemble transform kalman filter. part i: Theoretical aspects.  Monthly Weather Review, 129(3):420–436, Mar 2001.
  * [3] M. Bocquet.  Ensemble kalman filtering without the intrinsic need for inflation.  Nonlin. Processes Geophys., 18:735–750, 2011.
  * [4] S. A. Bourdarie and V. F. Maget.  Electron radiation belt data assimilation with an ensemble kalman filter relying on the salammbô code.  Annales Geophysicae, 30(6):929–943, jun 2012.
  * [5] F. Bouttier.  A dynamical estimation of error covariances in an assimilation system.  Monthly Weather Review, 122:2376–2390, 1994.
  * [6] M. Buehner.  Ensemble-derived stationary and flow-dependent background-error covariances: Evaluation in a quasi-operational nwp setting.  Quarterly Journal of the Royal Meteorological Society, 131:1013–1043, 2005.
  * [7] G. Burgers, P. J. van Leeuwen, and G. Evensen.  Analysis scheme in the ensemble kalman filter.  Monthly Weather Review, 126:1719–1724, 1998.
  * [8] CAMS.  The cams european air quality ensemble forecasts welcomes two new state-of-the-art models, 2022.
  * [9] Stephen E. Cohn.  Dynamics of short-term univariate forecast error covariances.  Monthly Weather Review, 121:3123–3149, 1993.
  * [10] Nour Dahmen, François Rogier, and Vincent Maget.  On the modelling of highly anisotropic diffusion for electron radiation belt dynamic codes.  Computer Physics Communications, 254:107342, sep 2020.
  * [11] J Derber, R. Purser, W.-S Wu, R Treadon, Manuel Pondeca, D Parrish, and Daryl Kleist.  Flow dependent jb in a global grid-point 3d-var.  09 2003.
  * [12] David C. Dowell, Fuqing Zhang, Louis J. Wicker, Chris Snyder, and N. Andrew Crook.  Wind and temperature retrievals in the 17 may 1981 arcadia, oklahoma, supercell: Ensemble kalman filter experiments.  Monthly Weather Review, 132(8):1982–2005, aug 2004.
  * [13] Ronald M. Errico, Tomislava Vukićević, and Kevin Raeder.  Comparison of initial and lateral boundary condition sensitivity for a limited-area model.  Tellus A: Dynamic Meteorology and Oceanography, 45(5):539, jan 1993\.
  * [14] G. Evensen.  Data Assimilation: The Ensemble Kalman Filter.  Springer-Verlag Berlin Heidelberg, 2009.
  * [15] Shay Gilpin, Tomoko Matsuo, and Stephen E. Cohn.  Continuum covariance propagation for understanding variance loss in advective systems.  SIAM/ASA Journal on Uncertainty Quantification, 10(3):886–914, 2022\.
  * [16] Nils Gustafsson.  Control of lateral boundary conditions in four-dimensional variational data assimilation for a limited area model.  Tellus A: Dynamic Meteorology and Oceanography, 64(1):17518, mar 2012.
  * [17] J. Harlim and B. R. Hunt.  A non-gaussian ensemble filter for assimilating infrequent noisy observations.  Tellus A, 59:225–237, 2007.
  * [18] P. Houtekamer and H. Mitchell.  A sequential ensemble kalman filter for atmospheric data assimilation.  Monthly Weather Review, 129:123–137, 2001.
  * [19] P. Houtekamer and F. Zhang.  Review of the ensemble kalman filter for atmospheric data assimilation.  Monthly Weather Review, 144:4489–4532, 2016.
  * [20] A. Jazwinski.  Stochastic Processes and Filtering Theory.  Academic Press, New York, 1970.
  * [21] A. Lorenc.  The potential of ensemble kalman filter for nwp – a comparison with 4d-var.  Quarterly Journal of the Royal Meteorological Society, 129:3183–3203, 2003.
  * [22] V. Maget, A. Sicard-Piet, S. Bourdarie, D. Lazaro, D. L. Turner, I. A. Daglis, and I. Sandberg.  Improved outer boundary conditions for outer radiation belt data assimilation using THEMIS-SST data and the Salammbo-EnKF code.  Journal of Geophysical Research: Space Physics, 120(7):2015JA021001, July 2015.
  * [23] V. Marécal, V.-H. Peuch, C. Andersson, S. Andersson, J. Arteta, M. Beekmann, A. Benedictow, R. Bergström, B. Bessagnet, A. Cansado, F. Chéroux, A. Colette, A. Coman, R. L. Curier, H. A. C. Denier van der Gon, A. Drouin, H. Elbern, E. Emili, R. J. Engelen, H. J. Eskes, G. Foret, E. Friese, M. Gauss, C. Giannaros, J. Guth, M. Joly, E. Jaumouillé, B. Josse, N. Kadygrov, J. W. Kaiser, K. Krajsek, J. Kuenen, U. Kumar, N. Liora, E. Lopez, L. Malherbe, I. Martinez, D. Melas, F. Meleux, L. Menut, P. Moinat, T. Morales, J. Parmentier, A. Piacentini, M. Plu, A. Poupkou, S. Queguiner, L. Robertson, L. Rouïl, M. Schaap, A. Segers, M. Sofiev, L. Tarasson, M. Thomas, R. Timmermans, Á. Valdebenito, P. van Velthoven, R. van Versendaal, J. Vira, and A. Ung.  A regional air quality forecasting system over europe: the MACC-II daily ensemble production.  Geoscientific Model Development, 8(9):2777–2813, sep 2015.
  * [24] S. Massart, B. Pajot, A. Piacentini, and O. Pannekoucke.  On the merits of using a 3d-fgat assimilation scheme with an outer loop for atmospheric situations governed by transport.  Monthly Weather Review, 138:4509–4522, 2010.
  * [25] Richard Ménard, Sergey Skachko, and Olivier Pannekoucke.  Numerical discretization causing error variance loss and the need for inflation.  Quarterly Journal of the Royal Meteorological Society, aug 2021\.
  * [26] I. Mirouze and A. T. Weaver.  Representation of correlation functions in variational assimilation using an implicit diffusion operator.  Quarterly Journal of the Royal Meteorological Society, 136(651):1421–1443, Jul 2010.
  * [27] P. Nutter, D. Stensrud, and M. Xue.  Effects of coarsely resolved and temporally interpolated lateral boundary conditions on the dispersion of limited-area ensemble forecasts.  Monthly Weather Review, 132(10):2358–2377, Oct 2004.
  * [28] P. Nutter, M. Xue, and D. Stensrud.  Application of lateral boundary condition perturbations to help restore dispersion in limited-area ensemble forecasts.  Monthly Weather Review, 132(10):2378–2390, Oct 2004.
  * [29] Christopher Paciorek and Mark Schervish.  Spatial modelling using a new class of nonstationary covariance functions.  Environmetrics, 17:483–506, 08 2006.
  * [30] O Pannekoucke.  opannekoucke/sympkf: A symbolic and computational toolbox for the design of parametric kalman filter dynamics, 2021.
  * [31] O. Pannekoucke, L. Berre, and G. Desroziers.  Background error correlation length-scale estimates and their sampling statistics.  Quaterly Journal Royal Meteorological Society, 134:497–508, 2008\.
  * [32] O. Pannekoucke, M. Bocquet, and R. Ménard.  Parametric covariance dynamics for the nonlinear diffusive burgers’ equation.  Nonlinear Processes in Geophysics, 2018:1–21, 2018.
  * [33] O. Pannekoucke and S. Massart.  Estimation of the local diffusion tensor and normalization for heterogeneous correlation modelling using a diffusion equation.  Quarterly Journal of the Royal Meteorological Society., 134:1425–1438, 2008.
  * [34] O. Pannekoucke, R. Ménard, M. El Aabaribaoune, and M. Plu.  A methodology to obtain model-error covariances due to the discretization scheme from the parametric kalman filter perspective.  Nonlinear Processes in Geophysics, 28(1):1–22, 2021.
  * [35] Olivier Pannekoucke.  An anisotropic formulation of the parametric kalman filter assimilation.  Tellus A: Dynamic Meteorology and Oceanography, 73(1):1–27, 2021\.
  * [36] Olivier Pannekoucke and Philippe Arbogast.  SymPKF (v1.0): a symbolic and computational toolbox for the design of parametric kalman filter dynamics.  Geoscientific Model Development, 14(10):5957–5976, oct 2021.
  * [37] Olivier Pannekoucke and Ronan Fablet.  Pde-netgen 1.0: from symbolic partial differential equation (pde) representations of physical processes to trainable neural network representations.  Geoscientific Model Development, 13:3373–3382, 07 2020.
  * [38] Olivier Pannekoucke, Sophie Ricci, Sebastien Barthelemy, Richard Ménard, and Olivier Thual.  Parametric kalman filter for chemical transport models.  Tellus A: Dynamic Meteorology and Oceanography, 68(1):31547, 2016\.
  * [39] Antoine Perrot, Olivier Pannekoucke, and Vincent Guidard.  Toward a multivariate formulation of the PKF assimilation: application to a simplified chemical transport model (accepted in npg).  sep 2022.
  * [40] M. Raissi, P. Perdikaris, and G.E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, feb 2019.
  * [41] M. Sabathier, O Pannekoucke, and V Maget.  opannekoucke/pkf-boundary: Boundary conditions for the parametric kalman filter forecast, oct 2022.
  * [42] F. Scarpa and G. Milano.  Kalman smoothing technique applied to the inverse heat conduction problem.  Numerical Heat Transfer, Part B: Fundamentals, 28(1):79–96, jul 1995.
  * [43] Chris Snyder and Fuqing Zhang.  Assimilation of simulated doppler radar observations with an ensemble kalman filter.  Monthly Weather Review, 131(8):1663–1677, aug 2003.
  * [44] Volker Strassen.  Gaussian elimination is not optimal.  Numerische Mathematik, 13(4):354–356, aug 1969.
  * [45] R. Stull.  An introduction to boundary layer meteorology.  Kluwer Academic, 1988.
  * [46] M. Tippett, J. Anderson, C. Bishop, T. Hamill, and J. Whitaker.  Ensemble square root filters.  Monthly Weather Review, 131(7):1485–1490, Jul 2003.
  * [47] R. D. Torn, G. J. Hakim, and C. Snyder.  Boundary conditions for limited-area ensemble kalman filters.  Monthly Weather Review, 134:2490–2502, 2006.
  * [48] Sina Voshtani, Richard Ménard, Thomas W. Walker, and Amir Hakami.  Assimilation of GOSAT methane in the hemispheric CMAQ; part i: Design of the assimilation system.  Remote Sensing, 14(2):371, jan 2022.
  * [49] Tomislava Vukicevic and Ronald M. Errico.  The influence of artificial and physical factors upon predictability estimates using a complex limited-area model.  Monthly Weather Review, 118(7):1460–1482, jul 1990.
  * [50] Tomislava Vukicevic and Jan Paegle.  The influence of one-way interacting lateral boundary conditions upon predictability of flow in bounded numerical models.  Monthly Weather Review, 117(2):340–350, feb 1989.
  * [51] Thomas T. Warner, Lawrence E. Key, and Annette M. Lario.  Sensitivity of mesoscale-model forecast skill to some initial-data characteristics, data density, data position, analysis procedure and measurement error.  Monthly Weather Review, 117(6):1281–1310, jun 1989.
  * [52] A. Weaver and P. Courtier.  Correlation modelling on the sphere using a generalized diffusion equation (tech. memo. ecmwf, num. 306).  Quarterly Journal of the Royal Meteorological Society, 127:1815–1846, 2001.