Statistical solutions of hyperbolic systems of conservation laws: numerical approximation

TL;DR

This paper introduces a numerical algorithm combining finite volume methods and Monte Carlo sampling to approximate statistical solutions of hyperbolic conservation laws, with proven convergence and demonstrated properties.

Contribution

It presents a new computational approach for approximating statistical solutions, integrating high-resolution finite volume methods with Monte Carlo sampling, and proves convergence under certain conditions.

Findings

Algorithm converges to statistical solutions

Numerical experiments confirm theoretical convergence

Reveals properties of statistical solutions

Abstract

Statistical solutions are time-parameterized probability measures on spaces of integrable functions, that have been proposed recently as a framework for global solutions and uncertainty quantification for multi-dimensional hyperbolic system of conservation laws. By combining high-resolution finite volume methods with a Monte Carlo sampling procedure, we present a numerical algorithm to approximate statistical solutions. Under verifiable assumptions on the finite volume method, we prove that the approximations, generated by the proposed algorithm, converge in an appropriate topology to a statistical solution. Numerical experiments illustrating the convergence theory and revealing interesting properties of statistical solutions, are also presented.

Figures

Equations

$${\partial }_{t}u+{\nabla }_x\cdot f(u)$$

$$u(x,0)$$

$${\int }_{{\mathbb{R}}_{+}}{\int }_{{\mathbb{R}}^d}{\partial }_t\phi (x,t)u(x,t)+{\nabla }_x\phi (x,t)\cdot f(u(x,t))dxdt+{\int }_{{\mathbb{R}}^d}\phi (x,0)\overline{u}(x)dx=0$$

$${\int }_{{\mathbb{R}}_{+}}{\int }_{{\mathbb{R}}^d}{\partial }_t\phi (x,t)\eta (u(x,t))+{\nabla }_x\phi (x,t)\cdot q(u(x,t))dxdt+{\int }_{{\mathbb{R}}^d}\phi (x,0)\eta (\bar{u}(x))dx⩾0$$

$${\int }_{{\mathbb{R}}^d}|u(x,t){|}^{p}dx⩽{\int }_{{\mathbb{R}}^d}|\bar{u}(x){|}^{p}dx\forall t>0,$$

$$\bigl{\langle}\mu_{n},F\bigr{\rangle}\to\bigl{\langle}\mu,F\bigr{\rangle}\qquad\text{as }n\to\infty$$

$$L_g(u):={\int }_{D^k}g(x_1,\dots ,x_k,u(x_1),\dots ,u(x_k))dx.$$

$$|g(x,\xi )|⩽\sum _{\alpha \in \{0,1{\}}^k}{\phi }_{|\bar{\alpha }|}(x_{\bar{\alpha }})|{\xi }^{\alpha }{|}^p\forall x\in D^k,\xi \in U^k$$

$$\begin{split}\big{|}g(x,\zeta)-g(y,\xi)\big{|}\leqslant\sum_{i=1}^{k}|\zeta_{i}-\xi_{i}|\max\big{(}|\xi_{i}|,|\zeta_{i}|\big{)}^{p-1}h(\hat{x}^{i},\hat{\xi}^{i})\end{split}$$

$$\displaystyle\mathcal{C}^{p}(D;U):=\big{\{}L_{g}\ :\ g\in\mathcal{H}^{k,p}(D;U),\ k\in\mathbb{N}\big{\}}$$

$$\displaystyle\mathcal{C}^{p}_{1}(D;U):=\big{\{}L_{g}\ :\ g\in\mathcal{H}^{k,p}_{1}(D;U),\ k\in\mathbb{N}\big{\}}$$

$$|g(x,\xi )|⩽{\phi }_1(x)+{\phi }_0|\xi {|}^p$$

$$|g(x_1,x_2,{\xi }_1,{\xi }_2)|⩽{\phi }_2(x_1,x_2)+{\phi }_1(x_1)|{\xi }_2{|}^p+{\phi }_1(x_2)|{\xi }_1{|}^p+{\phi }_0|{\xi }_1{|}^p|{\xi }_2{|}^p$$

$$\bigl{\langle}\nu^{k},g\bigr{\rangle}_{\mathcal{H}^{k}}=\int_{D^{k}}\bigl{\langle}\nu^{k}_{x},g(x)\bigr{\rangle}\,dx$$

$$\int_{D}\bigl{\langle}\nu^{1}_{x},|\xi|^{p}\bigr{\rangle}\,dx<+\infty.$$

$$\omega_{r}^{p}(\nu^{2}):=\int_{D}\mathchoice{{\vbox{\hbox{$\textstyle-$ }}\kern-7.31105pt}}{{\vbox{\hbox{$\scriptstyle-$ }}\kern-5.50664pt}}{{\vbox{\hbox{$\scriptscriptstyle-$ }}\kern-4.4333pt}}{{\vbox{\hbox{$\scriptscriptstyle-$ }}\kern-3.96664pt}}\!\int_{B_{r}(x)}\bigl{\langle}\nu^{2}_{x,y},|\xi_{1}-\xi_{2}|^{p}\bigr{\rangle}\,dydx.$$

$${\int }_{L^p}\Vert u{\Vert }_{L^p}^{p}d\mu (u)<\infty$$

$$\int_{D^{k}}\bigl{\langle}\nu^{k},g(x)\bigr{\rangle}\,dx=\int_{L^{p}}\int_{D^{k}}g(x,u(x))\,dxd\mu(u)$$

$$\operatorname*{ess\,sup}_{u\in L^{p}}\|u\|_{L^{p}}=\limsup_{k\to\infty}\left(\int_{D^{k}}\bigl{\langle}\nu^{k}_{x},|\xi_{1}|^{p}\cdots|\xi_{k}|^{p}\bigr{\rangle}\,dx\right)^{1/kp}$$

$$\underset{u\in L^p}{\mathrm{esssup}}\Vert u{\Vert }_{L^p(D;U)}^p$$

$$=\underset{k\to \infty }{\lim }{\left({\int }_{L^p(D;U)}{\int }_{D^k}|u(x_1){|}^p\cdots |u(x_k){|}^{p}dxd\mu (u)\right)}^{1/k}$$

$$\displaystyle=\lim_{k\to\infty}\left(\int_{D^{k}}\bigl{\langle}\nu^{k}_{x},|\xi_{1}|^{p}\cdots|\xi_{k}|^{p}\bigr{\rangle}\,dx\right)^{1/k}.$$

$$\limsup_{k\to\infty}\left(\int_{D^{k}}\bigl{\langle}\nu^{k}_{x},|\xi_{1}|^{p}\cdots|\xi_{k}|^{p}\bigr{\rangle}\,dx\right)^{1/kp}\leqslant M.$$

$$\displaystyle\sup_{n\in\mathbb{N}}\bigl{\langle}\nu^{1}_{n},|\xi|^{p}\bigr{\rangle}_{\mathcal{H}^{1}}$$

$$\displaystyle\lim_{r\to 0}\limsup_{n\to\infty}\omega_{r}^{p}\bigl{(}\nu^{2}_{n}\bigr{)}$$

$$\bigl{\langle}\nu^{k},g\bigr{\rangle}_{\mathcal{H}^{k}}\leqslant\liminf_{j\to\infty}\bigl{\langle}\nu_{n_{j}}^{k},g\bigr{\rangle}_{\mathcal{H}^{k}}.$$

$$\Vert u{\Vert }_{L^p}⩽M\text{for μn-a.e. }u\in L^p(D;U)\text{ for every }n\in \mathbb{N}.$$

$$\lim_{j\to\infty}\int_{D^{k}}\big{|}\bigl{\langle}\nu^{k}_{n_{j},x},g(x)\bigr{\rangle}-\bigl{\langle}\nu^{k}_{x},g(x)\bigr{\rangle}\big{|}\,dx=0$$

$$|g(t,x,\xi )|⩽\sum _{\alpha \in \{0,1{\}}^k}{\phi }_{|\bar{\alpha }|}(t,x_{\bar{\alpha }})|{\xi }^{\alpha }{|}^p\forall x\in D^k,\xi \in U^k$$

$$\begin{split}\big{|}g(t,x,\zeta)-g(t,y,\xi)\big{|}\leqslant\psi(t)\sum_{i=1}^{k}|\zeta_{i}-\xi_{i}|\max\big{(}|\xi_{i}|,|\zeta_{i}|\big{)}^{p-1}h(t,\hat{x}^{i},\hat{\xi}^{i})\end{split}$$

$$\bigl{\langle}\nu^{k},g\bigr{\rangle}_{\mathcal{H}^{k}}=\int_{0}^{T}\bigl{\langle}\rho(t),g(t,\cdot)\bigr{\rangle}_{\mathcal{H}^{k}}\,dt\qquad\forall\ g\in\mathcal{H}^{k}_{0}([0,T),D;U).$$

$$\operatorname*{ess\,sup}_{t\in[0,T)}\int_{D}\bigl{\langle}\nu^{1}_{t,x},|\xi|^{p}\bigr{\rangle}\,dx\leqslant c^{p}<+\infty$$

$$\int_{0}^{T^{\prime}}\omega_{r}^{p}\big{(}\nu^{2}_{t}\big{)}\,dt\to 0\qquad\text{ as }r\to 0\text{ for all }T^{\prime}\in(0,T)$$

Full text

1. Introduction

Systems of conservation laws are a large class of nonlinear partial differential equations of the generic form

[TABLE]

Here, the unknown $u=u(x,t):D\times\mathbb{R}_{+}\to U$ is the vector of conserved variables and $f=(f^{1},\dots,f^{d}):\mathbb{R}^{N}\to\mathbb{R}^{N\times d}$ is the flux function. Here, we denote $\mathbb{R}_{+}:=[0,\infty)$ and $U:=\mathbb{R}^{N}$, and we let the physical domain $D\subset\mathbb{R}^{d}$ be some open, connected set. The system (1.1a) is hyperbolic if the flux Jacobian $\partial_{u}(f\cdot n)$ has real eigenvalues for all $n\in\mathbb{R}^{d}$ with $|n|=1$.

Many important models in physics and engineering are described by hyperbolic systems of conservation laws. Examples include the compressible Euler equations of gas dynamics, the shallow water equations of oceanography, the Magneto-Hydro-Dynamics (MHD) equations of plasma physics, and the equations of nonlinear elastodynamics [10].

1.1. Entropy Solutions.

It is well known that even if the initial data $\overline{u}$ in (1.1) is smooth, solutions of (1.1) develop discontinuities, such as shock waves and contact discontinuities, in finite time. Therefore, solutions to (1.1) are sought in the sense of distributions: A function $u\in L^{\infty}(\mathbb{R}^{d}\times\mathbb{R}_{+},\mathbb{R}^{N})$ is an weak solution of (1.1) if it satisfies

[TABLE]

for all test functions $\varphi\in C^{1}_{c}(\mathbb{R}^{d}\times\mathbb{R}_{+})$.

As weak solutions are not unique [10], it is necessary to augment them with additional admissibility criteria or entropy conditions to recover uniqueness. These entropy conditions are based on the existence of a so-called entropy pair — a pair of functions $\eta:\mathbb{R}^{N}\to\mathbb{R}$, $q:\mathbb{R}^{N}\to\mathbb{R}^{d}$, with $\eta$ convex and $q$ satisfying the compatibility condition $q^{\prime}=\eta^{\prime}\cdot f^{\prime}$ (where $f^{\prime}$ and $q^{\prime}$ are the Jacobian matrices of $f$ and $q$). An entropy solution of (1.1) is a weak solution that also satisfies the so-called entropy inequality

[TABLE]

for all nonnegative test functions $\varphi\in C^{1}_{c}(\mathbb{R}^{d}\times\mathbb{R}_{+})$. Depending on the availability of entropy pairs $(\eta,q)$, the entropy condition leads to various a priori bounds on $u$: If, say, $\eta(u)=|u|^{p}$ (or some perturbation thereof) for some $p\geqslant 1$ then (1.3) leads to

[TABLE]

see e.g. [10, 26, 29].

The global well-posedness of entropy solutions of (1.1) has been addressed both for (multi-dimensional) scalar conservation laws [32] and for systems in one space dimension (see [25, 4, 7, 29] and references therein). However, there are no global existence results for entropy solutions of multi-dimensional systems of conservation laws with generic initial data. On the other hand, it has been established recently in [11, 9] that entropy solutions for some systems of conservation laws (such as isentropic Euler equations in two space dimensions) may not be unique. This is a strong indication that the paradigm of entropy solutions is not the correct framework for the well-posedness of multi-dimensional systems of hyperbolic conservation laws.

1.2. Numerical schemes

A wide variety of numerical methods have been developed to approximate entropy solutions of (1.1) in a robust and efficient manner. These include finite volume, (conservative) finite difference, discontinuous Galerkin (DG) finite element and spectral (viscosity) methods, see textbooks [26, 28] for further details. Rigorous convergence results of numerical methods to entropy solutions are only available for scalar conservation laws (see e.g. [26] for monotone schemes and [14] for high-order schemes) and for some specific numerical methods for one-dimensional systems ([25] for Glimm’s scheme and [29] for front tracking).

There are no rigorous convergence results to entropy solutions for any numerical schemes approximating multi-dimensional systems of conservation laws. To the contrary, several numerical experiments, such as those presented recently in [15, 19], strongly suggest that there is no convergence of approximations generated by standard numerical schemes for (1.1), as the mesh is refined. This has been attributed to the emergence of turbulence-like structures at smaller and smaller scales upon mesh refinement (see Figure 4 of [15]).

1.3. Measure-valued and Statistical solutions.

Given the lack of well-posedness of entropy solutions for multi-dimensional systems of conservation laws and the lack of convergence of numerical approximations to them, it is natural to seek alternative solution paradigms for (1.1). A possible solution framework is that of entropy measure-valued solutions, first proposed by DiPerna in [12], see also [13]. Measure-valued solutions are Young measures [48], that is, space-time parameterized probability measures on phase space $\mathbb{R}^{N}$ of (1.1). Global existence of entropy measure-valued solutions has been considered in [12, 8] and in [15, 19] where the authors constructed entropy measure-valued solutions by proving convergence of a Monte Carlo type ensemble-averaging algorithm, based on underlying entropy stable finite difference schemes.

Although entropy measure-valued solutions for multi-dimensional systems of conservation laws exist globally, it is well known that they are not necessarily unique; see [42, 19] and references therein. In particular, one can even construct multiple entropy measure-valued solutions for scalar conservation laws for the same measure-valued initial data [19]. Although generic measure-valued solutions might not be unique, numerical experiments presented in [15] indicate that measure-valued solutions of (1.1), computed with the ensemble-averaging algorithm of [15], are stable with respect to initial perturbations and to the choice of underlying numerical method. This suggests imposing additional constraints on entropy-measure valued solutions in order to recover uniqueness.

In [16], the authors implicated the lack of information about (multi-point) statistical correlations in Young measures as a possible cause of the non-uniqueness of entropy measure-valued solutions. Consequently, they introduced a stronger solution paradigm termed statistical solutions for hyperbolic systems of conservation laws (1.1). Statistical solutions are time-parameterized probability measures on some Lebesgue space $L^{p}(D;U)$ satisfying (1.1a) in an averaged sense. The choice of the exponent $p\geqslant 1$ depends on the available a priori bounds for solution of (1.1), such as (1.4). It was shown in [16] that probability measures on $L^{p}(D;U)$ can be identified with (and indeed are equivalent to) correlation measures — a hierarchy of Young measures defined on tensorized versions of the domain $D$ and the phase space $U$ in (1.1). Statistical solutions have also been introduced in the context of the incompressible Navier–Stokes equations by Foiaş et al.; see [21] and references therein.

In [16], the authors defined statistical solutions of systems of conservation laws (1.1a) by requiring that the moments of the time-parameterized probability measure on $L^{p}(D)$ (or equivalently, of the underlying correlation measure) satisfy an infinite set of (tensorized) partial differential equations (PDEs), consistent with (1.1a).

The first member of the hierarchy of correlation measures for a statistical solution, is a (classical) Young measure and it can be shown to be an entropy measure valued solution of (1.1), in the sense of DiPerna [12]. The $k$-th member ($k\geqslant 2$) of the hierarchy represents $k$-point spatial correlations. Thus, a statistical solution can be thought of as a measure-valued solution, augmented with information about all possible (multi-point) spatial correlations [16]. Consequently, statistical solutions contain much more information than measure-valued solutions.

In [16], the authors constructed a canonical statistical solution for scalar conservation laws in terms of the data-to-solution semi-group of Kruzkhov [32] and showed that this statistical solution is unique under a suitable entropy condition. Numerical approximation of statistical solutions of scalar conservation laws was considered in [17], where the authors proposed Monte Carlo and multi-level Monte Carlo (MLMC) algorithms to compute statistical solutions and showed their convergence.

1.4. Aims and scope of this paper.

Given this background, our main aim in this paper is to study statistical solutions for multi-dimensional systems of conservation laws. To this end, we obtain the following results:

• •

We propose a Monte Carlo ensemble averaging based algorithm for computing statistical solutions of systems of conservation laws. This algorithm is a variant of the Monte Carlo algorithms presented in [15] and [17].

• •

Under reasonable assumptions on the underlying numerical scheme, we prove convergence of the ensemble-averaging algorithm to a statistical solution. It is highly non-trivial to identify an appropriate topology on time parameterized probability measures on $L^{p}(D)$ in order to prove convergence of the computed statistical solutions. To this end, we find a suitable topology and prescribe novel sufficient conditions that ensure convergence in this topology.

• •

We present several numerical experiments that illustrate the robustness of our proposed algorithm and also reveal interesting properties of statistical solutions of (1.1a).

As a consequence of our convergence theorem, we establish a conditional global existence result for multi-dimensional systems of conservation laws. Moreover, we also propose an entropy condition under which we prove a weak-strong uniqueness result for statistical solutions, that is, we prove that if there exists a statistical solution of sufficient regularity (in a sense made precise in Section 3), then all entropy statistical solutions agree with it.

The rest of the paper is organized as follows: In Section 2, we provide the mathematical framework by describing the concepts of correlation measures and statistical solutions. We also provide characterizations of the topology on probability measures on $L^{p}(D)$, in which our subsequent numerical approximations will converge. The entropy condition and the weak-strong uniqueness of statistical solutions are presented in Section 3 and the Monte Carlo ensemble-averaging algorithm (and its convergence) is presented in Section 4. Numerical experiments are presented in Section 5 and the results of the paper are summarized and discussed in Section 6.

2. Probability measures on $L^{p}(D;U)$ and Statistical solutions

In the usual, deterministic interpretation of (1.1a), one attempts to find a function $u=u(t):D\to U$ satisfying (1.1a) in a weak or strong sense. (Here, as in the Introduction, we let $D\subset\mathbb{R}^{d}$ be an open, connected set and we denote $U:=\mathbb{R}^{N}$.) By contrast, a statistical solution of (1.1a) is a probability measure $\mu=\mu_{t}$ distributed over such functions $u$ and satisfying (1.1a) in some averaged sense. Solutions of (1.1a) are most naturally found in (some subspace of) $L^{p}(D;U)$, so $\mu_{t}$ is required to be a probability measure on $L^{p}(D;U)$ at each time $t$. In order to write down constitutive equations for $\mu$, it is more natural to work with finite-dimensional projections or marginals of $\mu$; these are the so-called correlation measures [16]. In this section we provide a self-contained description of correlation measures, probability measures over $L^{p}(D;U)$, and statistical solutions of (1.1a).

In order to link probability measures on $L^{p}$ to their finite-dimensional marginals, we prove in Section 2.1 that a sequence of such measures converges weakly if and only if it converges with respect to a certain class $\mathcal{C}_{p}$ of finite-dimensional observables. In Section 2.2 we introduce correlation measures and we show that these are in a one-to-one relationship with probability measures over $L^{p}$, and that they are linked precisely through the finite-dimensional observables $\mathcal{C}_{p}$. We also prove a compactness result for sequences of correlation measures. In Section 2.3 we treat time-parametrized probability and correlation measures, and we prove measurability and compactness results. Finally, in Section 2.4 we provide the definition of statistical solutions of (1.1a).

For the sake of clarity, many of the proofs in this rather technical section have been moved to Appendices A, B and C.

Notation 2.1**.**

If $X$ is a topological space then we let $\mathcal{B}(X)$ denote the Borel $\sigma$-algebra on $X$, we let $\mathcal{M}(X)$ denote the set of signed Radon measures on $(X,\mathcal{B}(X))$, and we let $\mathcal{P}(X)\subset\mathcal{M}(X)$ denote the set of all probability measures on $(X,\mathcal{B}(X))$, i.e. all non-negative $\mu\in\mathcal{M}(X)$ with $\mu(X)=1$ (see e.g. [2, 5, 31]). For $k\in\mathbb{N}$ and a multiindex $\alpha\in\{0,1\}^{k}$ we write $|\alpha|=\alpha_{1}+\dots+\alpha_{k}$ and $\bar{\alpha}=\mathbbm{1}-\alpha=(1-\alpha_{1},\dots,1-\alpha_{k})$, and we let $x_{\alpha}$ be the vector of length $|\alpha|$ consisting of the elements $x_{i}$ of $x$ for which $\alpha_{i}$ is non-zero. For a vector $x=(x_{1},\dots,x_{k})$ we write $\hat{x}^{i}=(x_{1},\dots,x_{i-1},x_{i+1},\dots,x_{k})$. For a vector $\xi=(\xi_{1},\dots,\xi_{k})$ we write $|\xi^{\alpha}|=|\xi_{\alpha_{1}}|\cdots|\xi_{\alpha_{k}}|$ with the convention $\xi_{\alpha_{i}}=1$ if $\alpha_{i}=0$.

2.1. Probability measures on $L^{p}(D)$ and weak convergence

If $X$ is any topological space and we are given a sequence $\mu_{1},\mu_{2},\dots\in\mathcal{P}(X)$ and some $\mu\in\mathcal{P}(X)$, then we say that $\{\mu_{n}\}_{n\in\mathbb{N}}$ converges weakly to $\mu$, written $\mu_{n}\rightharpoonup\mu$, if

[TABLE]

for every $F\in C_{b}(X)$. (Here and elsewhere, \bigl{\langle}\mu,F\bigr{\rangle}=\int_{X}F(x)\,d\mu(x) denotes the expectation of $F$ with respect to $\mu$.) We will be particularly interested in the case $X=L^{p}(D;U)$, so to study weak convergence in this space we need to work with the space $C_{b}(L^{p}(D;U))$. In this section we will see that it is sufficient to prove (2.1) for a much smaller class of functionals $F$, namely those which depend only on finite-dimensional projections of $u\in L^{p}(D;U)$.

If $E$ and $V$ are Euclidean spaces then a measurable function $g:E\times V\to\mathbb{R}$ is called a Carathéodory function if $\xi\mapsto g(x,\xi)$ is continuous for a.e. $x\in E$ and $x\mapsto g(x,\xi)$ is measurable for every $\xi\in V$ (see e.g. [1, Section 4.10]). For a number $k\in\mathbb{N}$ and a Carathéodory function $g=g(x,\xi):D^{k}\times U^{k}\to\mathbb{R}$ we define the functional $L_{g}:L^{p}(D;U)\to\mathbb{R}$ by

[TABLE]

(Here, $D^{k}$ denotes the product space $D^{k}=D\times\dots\times D$, and similarly for $U^{k}$.) The above integral is clearly not well-defined for every Carathéodory function $g$, so we restrict our attention to the following class.

Definition 2.2**.**

For every $k\in\mathbb{N}$, we let $\mathcal{H}^{k,p}(D;U)$ denote the space of Carathéodory functions $g:D^{k}\times U^{k}\to\mathbb{R}$ satisfying

[TABLE]

for nonnegative functions $\varphi_{i}\in L^{1}(D^{i})$, $i=0,1,\dots,k$ (with the convention that $L^{1}(D^{0})\cong\mathbb{R}$; see also Example 2.3). We let $\mathcal{H}^{k,p}_{1}(D;U)\subset\mathcal{H}^{k,p}(D;U)$ denote the subspace of functions $g$ which are locally Lipschitz continuous, in the sense that there is some $r>0$ and some nonnegative $h\in\mathcal{H}^{k-1,p}(D;U)$ such that

[TABLE]

for every $x\in D^{k}$, $y\in B_{r}(x)$ and $\xi,\zeta\in U^{k}$. Last, we denote

[TABLE]

where $L_{g}$ is defined in (2.2).

Example 2.3**.**

For $k=1$ the condition (2.3) asserts that

[TABLE]

for $0\leqslant\varphi_{1}\in L^{1}(D)$ and $\varphi_{0}\in[0,\infty)$, and for $k=2$ that

[TABLE]

for $0\leqslant\varphi_{2}\in L^{1}(D^{2})$, $0\leqslant\varphi_{1}\in L^{1}(D)$ and $\varphi_{0}\in[0,\infty)$.

We will simply denote $\mathcal{H}^{k,p}=\mathcal{H}^{k,p}(D;U)$, etc. when the domain and image $D,U$ are clear from the context.

Lemma 2.4**.**

Every functional $L_{g}\in\mathcal{C}^{p}$ is well-defined and finite on $L^{p}(D;U)$. Every functional $L_{g}\in\mathcal{C}^{p}_{1}$ is continuous and is Lipschitz continuous on bounded subsets of $L^{p}(D;U)$.

Theorem 2.5**.**

Let $\mu_{n},\mu\in\mathcal{P}(L^{p}(D;U))$ for $n\in\mathbb{N}$ satisfy $\operatorname{supp}\mu,\operatorname{supp}\mu_{n}\subset B$ for all $n\in\mathbb{N}$, for some bounded set $B\subset L^{p}(D;U)$. Then $\mu_{n}\rightharpoonup\mu$ if and only if \bigl{\langle}\mu_{n},F\bigr{\rangle}\to\bigl{\langle}\mu,F\bigr{\rangle} for all $F\in\mathcal{C}^{p}_{1}(D;U)$.

The proofs of the above results can be found in Appendix A. The “only if” part of Theorem 2.5 is trivial, since every $F\in\mathcal{C}^{p}_{1}$ belongs to $C_{b}$ when restricted to a bounded set; the converse relies on an approximation argument found in [2].

2.2. Correlation measures

In short, a correlation measure prescribes the joint distribution of some uncertain quantity $u$ at any finite collection of spatial points $x_{1},\dots,x_{k}$. Below we provide the rigorous definition of correlation measures and then state the result from [16] on the equivalence between correlation measures and probability measures on $L^{p}(D;U)$.

We denote \mathcal{H}^{k}_{0}(D;U):=L^{1}\big{(}D^{k},C_{0}(U^{k})\big{)}. By identifying the expressions $g(x)(\xi)$ and $g(x,\xi)$, we can view $\mathcal{H}^{k}_{0}(D;U)$ as a subspace of $\mathcal{H}^{k,p}(D;U)$ for any $p\geqslant 1$ (with the choice $\varphi_{0},\dots,\varphi_{k-1}\equiv 0$ and $\varphi_{k}(x)=\|g(x)\|_{C_{0}(U^{k})}$ in (2.3)).

Theorem 2.6**.**

The dual of $\mathcal{H}^{k}_{0}(D;U)$ is the space \mathcal{H}^{k*}_{0}(D;U):=L^{\infty}_{w}\big{(}D^{k},\mathcal{M}(U^{k})\big{)}, the space of bounded, weak measurable maps from $D^{k}$ to $\mathcal{M}(U^{k})$, under the duality pairing*

[TABLE]

(where \bigl{\langle}\nu^{k}_{x},g(x)\bigr{\rangle}=\int_{U^{k}}g(x,\xi)\,d\nu^{k}_{x}(\xi) is the usual duality pairing between Radon measures $\mathcal{M}(U^{k})$ and continuous functions $C_{b}(U^{k})$).

For more details and references for the above result, see [3].

Definition 2.7** (Fjordholm, Lanthaler, Mishra [16]).**

A correlation measure is a collection $\bm{\nu}=(\nu^{1},\nu^{2},\dots)$ of maps $\nu^{k}\in\mathcal{H}^{k*}_{0}(D;U)$ satisfying for all $k=1,2,\dots$

1. (i)

$\nu^{k}_{x}\in\mathcal{P}(U^{k})$ for a.e. $x\in D^{k}$, and the map x\mapsto\bigl{\langle}\nu^{k}_{x},f\bigr{\rangle} is measurable for every $f\in C_{b}(U^{k})$. (In other words, $\nu^{k}$ is a Young measure from $D^{k}$ to $U^{k}$.) 2. (ii)

Symmetry: if $\sigma$ is a permutation of $\{1,\dots,k\}$ and $f\in C_{0}(U^{k})$ then \bigl{\langle}\nu^{k}_{\sigma(x)},f(\sigma(\xi))\bigr{\rangle}=\bigl{\langle}\nu^{k}_{x},f(\xi)\bigr{\rangle} for a.e. $x\in D^{k}$. 3. (iii)

Consistency: If $f\in C_{b}(U^{k})$ is of the form $f(\xi_{1},\dots,\xi_{k})=g(\xi_{1},\dots,\xi_{k-1})$ for some $g\in C_{0}(U^{k-1})$, then \bigl{\langle}\nu^{k}_{x_{1},\dots,x_{k}},f\bigr{\rangle}=\bigl{\langle}\nu^{k-1}_{x_{1},\dots,x_{k-1}},g\bigr{\rangle} Lebesgue-a.e. $(x_{1},\dots,x_{k})\in D^{k}$. 4. (iv)

$L^{p}$* integrability:*

[TABLE] 5. (v)

Diagonal continuity: $\lim_{r\to 0}\omega^{p}_{r}(\nu^{2})=0$, where

[TABLE]

Each element $\nu^{k}$ will be called a correlation marginal. The functional $\omega_{r}^{p}$ is called the modulus of continuity of $\bm{\nu}$. We let $\mathcal{L}^{p}(D;U)$ denote the set of all correlation measures.

The next result shows that there is a duality relation between correlation marginals and the probability measures $\mu\in\mathcal{P}(L^{p})$ discussed in the previous section.

Theorem 2.8** (Fjordholm, Lanthaler, Mishra [16]).**

For every correlation measure $\bm{\nu}\in\mathcal{L}^{p}(D;U)$ there is a unique probability measure $\mu\in\mathcal{P}(L^{p}(D;U))$ whose $p$-th moment is finite,

[TABLE]

and such that $\mu$ is dual to $\bm{\nu}$: the identity

[TABLE]

holds for every $g\in\mathcal{H}^{k}_{0}(D;U)$ and all $k\in\mathbb{N}$. Conversely, for every $\mu\in\mathcal{P}(L^{p}(D;U))$ satisfying (2.7) there is a unique correlation measure $\bm{\nu}\in\mathcal{L}^{p}(D;U)$ that is dual to $\mu$.

Remark 2.9**.**

By using Lebesgue’s dominated convergence theorem, it is not hard to show that the identity (2.8) can be extended to all $g\in\mathcal{H}^{k,p}(D;U)$, as long as both integrals are well-defined. In particular, this is true if $\mu$ is supported on a bounded subset of $L^{p}(D;U)$.

Later on, we will be particularly interested in those $\mu\in\mathcal{P}(L^{p})$ that have bounded support. The following lemma shows how the property of having bounded support can be expressed in terms of the corresponding correlation measure.

Lemma 2.10**.**

Let $\bm{\nu}\in\mathcal{L}^{p}(D;U)$ and $\mu\in\mathcal{P}(L^{p}(D;U))$ be dual to one another. Then

[TABLE]

where the “$\operatorname*{ess\,sup}$” is taken with respect to $\mu$.

Proof.

From the identity $\|f\|_{L^{\infty}(X;\mu)}=\lim_{k\to\infty}\|f\|_{L^{k}(X;\mu)}$, valid for any finite measure $\mu$, we get

[TABLE]

∎

Definition 2.11**.**

We let $\mathcal{L}^{p}_{b}(D;U)$ denote the subset of correlation measures $\bm{\nu}\in\mathcal{L}^{p}(D;U)$ with bounded support, in the sense that there is an $M>0$ such that

[TABLE]

Definition 2.12**.**

If $\bm{\nu}_{n},\bm{\nu}\in\mathcal{L}^{p}(D;U)$ for $n\in\mathbb{N}$ then we say that $\bm{\nu}_{n}$ converges weak* to $\bm{\nu}$ as $n\to\infty$ (written $\bm{\nu}_{n}\overset{*}{\rightharpoonup}\bm{\nu}$) if $\nu^{k}_{n}\overset{*}{\rightharpoonup}\nu^{k}$ as $n\to\infty$, that is, if \bigl{\langle}\nu^{k}_{n},g\bigr{\rangle}_{\mathcal{H}^{k}}\to\bigl{\langle}\nu^{k},g\bigr{\rangle}_{\mathcal{H}^{k}} for all $g\in\mathcal{H}^{k}_{0}(D;U)$ and all $k\in\mathbb{N}$.

If $\bm{\nu}_{n},\bm{\nu}\in\mathcal{L}^{p}_{b}(D;U)$ for $n\in\mathbb{N}$ then we say that $(\bm{\nu}_{n})_{n\in\mathbb{N}}$ converges weakly to $\bm{\nu}$ as $n\to\infty$ (written $\bm{\nu}_{n}\rightharpoonup\bm{\nu}$) if \bigl{\langle}\nu^{k}_{n},g\bigr{\rangle}_{\mathcal{H}^{k}}\to\bigl{\langle}\nu^{k},g\bigr{\rangle}_{\mathcal{H}^{k}} for every $g\in\mathcal{H}^{k,p}_{1}(D;U)$.

Note that $\bm{\nu}\in\mathcal{L}^{p}_{b}$ implies that \bigl{\langle}\nu^{k},g\bigr{\rangle}_{\mathcal{H}^{k}} is well-defined and finite for any $g\in\mathcal{H}^{k,p}$ (cf. Definition 2.2).

We next show a compactness result which can be thought of as Kolmogorov’s compactness theorem (cf. [29, Theorem A.5]) for correlation measures.

Theorem 2.13**.**

Let $\bm{\nu}_{n}\in\mathcal{L}^{p}(D;U)$ for $n=1,2,\dots$ be a sequence of correlation measures such that

[TABLE]

for some $c>0$ (where $\omega_{r}^{p}$ is defined in Definition 2.7(v)). Then there exists a subsequence $(n_{j})_{j=1}^{\infty}$ and some $\bm{\nu}\in\mathcal{L}^{p}(D;U)$ such that

1. (i)

$\bm{\nu}_{n_{j}}\overset{*}{\rightharpoonup}\bm{\nu}$* as $j\to\infty$, that is, \bigl{\langle}\nu^{k}_{n_{j}},g\bigr{\rangle}_{\mathcal{H}^{k}}\to\bigl{\langle}\nu^{k},g\bigr{\rangle}_{\mathcal{H}^{k}} for every $g\in\mathcal{H}^{k}_{0}(D;U)$ and every $k\in\mathbb{N}$* 2. (ii)

\bigl{\langle}\nu^{1},|\xi|^{p}\bigr{\rangle}_{\mathcal{H}^{1}}\leqslant c^{p}** 3. (iii)

$\omega_{r}^{p}(\nu^{2})\leqslant\liminf_{n\to\infty}\omega_{r}^{p}(\nu^{2}_{n})$* for every $r>0$* 4. (iv)

for $k\in\mathbb{N}$, let $\varphi\in L^{1}_{{\text{loc}}}(D^{k})$ and $\kappa\in C(U^{k})$ be nonnegative, and let $g(x,\xi):=\varphi(x)\kappa(\xi)$. Then

[TABLE] 5. (v)

Assume moreover that the domain $D\subset\mathbb{R}^{d}$ is bounded and that $\bm{\nu}_{n}$ have uniformly bounded support, in the sense that (2.10) holds for all $\bm{\nu}_{n}$ for a fixed $M>0$, or equivalently,

[TABLE]

Then observables converge strongly:

[TABLE]

*for every $g\in\mathcal{H}^{k,p}_{1}(D;U)$. In particular, $\bm{\nu}_{n_{j}}\rightharpoonup\bm{\nu}$. *

The proof is given in Appendix B.

Remark 2.14**.**

(2.15) implies in particular that \bigl{\langle}\nu^{k}_{n_{j}},g\bigr{\rangle}_{\mathcal{H}^{k}}=\bigl{\langle}\mu_{n_{j}},L_{g}\bigr{\rangle} converges for any $g\in\mathcal{H}^{k,p}_{1}$, where $\mu_{n}\in\mathcal{P}(L^{p})$ is dual to $\bm{\nu}_{n}$ (see Theorem 2.8). By Theorem 2.5, this is equivalent to saying that $\mu_{n_{j}}$ converges weakly to $\mu$. Since, by hypothesis, the $p$th moment of $\mu_{n}$ is uniformly bounded, the sequence $\mu_{n}$ converges to $\mu$ in the Wasserstein distance; see Definition 3.1 and [45, Chapter 7].

Remark 2.15**.**

Theorem 2.13 can most likely be extended to provide a complete characterization of compact subsets of $\mathcal{L}^{p}(D;U)$. Since we only require sufficient conditions for compactness, we do not pursue this generalization here.

2.3. Time-parameterized probability measures on $L^{p}$

Let $T\in(0,\infty]$. To take into account the evolutionary nature of the PDE (1.1a), we will add time-dependence to the probability measures considered in Section 2.1 by considering maps $\mu:[0,T)\to\mathcal{P}(L^{p}(D;U))$. Note the distinction between time-parametrized maps $\mu:[0,T)\to\mathcal{P}(L^{p}(D;U))$ and probability measures $\gamma$ on, say, the space $L^{\infty}([0,T);L^{p}(D;U))$. Every such measure $\gamma$ would correspond to a unique $\mu$, but not vice versa; when “projecting” $\gamma$ onto $\mu$, any information about correlation between function values $u(t_{1})$, $u(t_{2})$ at different times $t_{1}$, $t_{2}$ is lost. Given the evolutionary nature of the PDE (1.1a), we have chosen to work with “$\mu$” measures in order to preserve the direction of time in the underlying PDE.

Notation 2.16**.**

We denote the set of Carathéodory functions depending on space and time by $\mathcal{H}^{k}_{0}([0,T),D;U):=L^{1}([0,T)\times D^{k};C_{0}(U^{k}))$ and its dual space by $\mathcal{H}^{k*}_{0}([0,T),D;U):=L^{\infty}_{w}([0,T)\times D^{k};\mathcal{M}(U^{k}))$.

Analogously to Definition 2.2, we let $\mathcal{H}^{k,p}([0,T),D;U)$ denote the space of Carathéodory functions $g:[0,T)\times D^{k}\times U^{k}\to\mathbb{R}$ satisfying

[TABLE]

for nonnegative functions $\varphi_{i}\in L^{\infty}([0,T);L^{1}(D^{i}))$, $i=0,1,\dots,k$. We let $\mathcal{H}^{k,p}_{1}([0,T),D;U)\subset\mathcal{H}^{k,p}([0,T),D;U)$ denote the subspace of functions $g$ satisfying the local Lipschitz condition

[TABLE]

for every $x\in D^{k}$, $y\in B_{r}(x)$ for some $r>0$, for some nonnegative $h\in\mathcal{H}^{k-1,p}([0,T),D;U)$ and $0\leqslant\psi(t)\in L^{\infty}([0,T))$.

The following lemma shows that it is meaningful to “evaluate” an element $\nu^{k}\in\mathcal{H}^{k*}_{0}([0,T),D;U)$ at (almost) any time $t\in[0,T)$.

Lemma 2.17**.**

Let $\nu^{k}\in\mathcal{H}^{k*}_{0}([0,T),D;U)$. Then there exists a map $\rho:[0,T)\to\mathcal{H}^{k*}_{0}(D;U)$, uniquely defined for a.e. $t\in[0,T)$, such that t\mapsto\bigl{\langle}\rho(t),g\bigr{\rangle}_{\mathcal{H}^{k}_{0}} is measurable for all $g\in\mathcal{H}^{k}_{0}(D;U)$, and

[TABLE]

The proof of this lemma is given in Appendix C. Henceforth we will not make distinctions between these two representations of elements of $\mathcal{H}^{k*}_{0}([0,T),D;U)$, and denote them both by $\nu^{k}$.

Definition 2.18**.**

A time-dependent correlation measure is a collection $\bm{\nu}=(\nu^{1},\nu^{2},\dots)$ of functions $\nu^{k}\in\mathcal{H}^{k*}_{0}([0,T),D;U)$ such that

1. (i)

$(\nu^{1}_{t},\nu^{2}_{t},\dots)\in\mathcal{L}^{p}(D;U)$ for a.e. $t\in[0,T)$ 2. (ii)

$L^{p}$* integrability:*

[TABLE] 3. (iii)

Diagonal continuity (DC):

[TABLE]

where $\omega_{r}^{p}$ was defined in (2.6).

We denote the set of all time-dependent correlation measures by $\mathcal{L}^{p}([0,T),D;U)$.

Remark 2.19**.**

By Lemma 2.17, the objects $\nu^{k}_{t}$ are well-defined for a.e. $t\in[0,T)$. Assertion (ii) requires that the $L^{p}$ bound should be uniform in $t$, and assertion (iii) requires that the modulus of continuity in the diagonal continuity requirement should be integrable in $t$.

Next, we prove a time-dependent version of the duality result Theorem 2.8.

Theorem 2.20**.**

For every time-dependent correlation measure $\bm{\nu}\in\mathcal{L}^{p}([0,T),D;U)$ there is a unique (up to subsets of $[0,T)$ of Lebesgue measure [math]) map $\mu:[0,T)\to\mathcal{P}(L^{p}(D;U))$ such that

1. (i)

the map

[TABLE]

is measurable for all $g\in\mathcal{H}^{k}_{0}(D;U)$, 2. (ii)

$\mu$* is $L^{p}$-bounded:*

[TABLE] 3. (iii)

$\mu$* is dual to $\bm{\nu}$: the identity*

[TABLE]

holds for a.e. $t\in[0,T)$, every $g\in\mathcal{H}^{k}_{0}(D;U)$ and all $k\in\mathbb{N}$.

Conversely, for every $\mu:[0,T)\to\mathcal{P}(L^{p}(D;U))$ satisfying (i) and (ii), there is a unique correlation measure $\bm{\nu}\in\mathcal{L}^{p}([0,T),D;U)$ satisfying (iii).

Proof.

Let $\bm{\nu}$ be given. Then for a.e. $t\in[0,T)$ we have $\bm{\nu}_{t}:=(\nu^{1}_{t},\nu^{2}_{t},\dots)\in\mathcal{L}^{p}(D;U)$, so by Theorem 2.8 there exists a unique $\mu_{t}\in\mathcal{P}(L^{p}(D;U))$ that is dual to $\bm{\nu}_{t}$, in the sense that (iii) holds. From the previous remark we know that t\mapsto\bigl{\langle}\nu^{k}_{t},g\bigr{\rangle}_{\mathcal{H}^{k}} is measurable for every $g\in\mathcal{H}^{k}_{0}(D;U)$, which (using (iii)) is precisely (i). Property (ii) follows by approximating $\xi\mapsto|\xi|^{p}$ by functions in $C_{0}(U)$.

Conversely, given $\mu$ satisfying (i) and (ii), Theorem 2.8 gives, for a.e. $t\in[0,T)$, the existence and uniqueness of $\bm{\nu}_{t}\in\mathcal{L}^{p}(D;U)$ satisfying (iii) as well as the $L^{p}$-bound (2.18). We claim that $(\bm{\nu}_{t})_{t\in[0,T)}$ defines a time-dependent correlation measure $\bm{\nu}\in\mathcal{L}^{p}([0,T),D;U)$. Indeed, define the linear functional $\nu^{k}$ by

[TABLE]

Then $\nu^{k}$ is well-defined on tensor product test functions $\theta(t)g(x)$, and

[TABLE]

Extending $\nu^{k}$ by linearity to all of $L^{1}([0,T)\times D^{k};C_{0}(U^{k}))$ produces a unique element $\nu^{k}\in L^{1}([0,T)\times D^{k};C_{0}(U^{k}))^{*}\cong L^{\infty}_{w}([0,T)\times D^{k};\mathcal{M}(U^{k}))$. Defining the collection $\bm{\nu}=(\nu^{1},\nu^{2},\dots)$, it only remains to show that $\nu^{2}$ satisfies the diagonal continuity requirement (2.19). Indeed, since

[TABLE]

for a.e. $t\in[0,T)$, the requirement (2.19) follows from the dominated convergence theorem. ∎

We denote the set of all maps $\mu:[0,T)\to\mathcal{P}(L^{p}(D;U))$ that are dual to some $\bm{\nu}\in\mathcal{L}^{p}([0,T),D;U)$ as $\mathcal{P}_{T}(L^{p}(D;U))$.

We conclude this section by proving a version of the compactness theorem for time-dependent correlation measures.

Theorem 2.21**.**

Let $\bm{\nu}_{n}\in\mathcal{L}^{p}([0,T),D;U)$ for $n=1,2,\dots$ be a sequence of correlation measures such that

[TABLE]

for some $c>0$ and all $T^{\prime}\in[0,T)$. Then there exists a subsequence $(n_{j})_{j=1}^{\infty}$ and some $\bm{\nu}\in\mathcal{L}^{p}([0,T),D;U)$ such that

1. (i)

$\bm{\nu}_{n_{j}}\overset{*}{\rightharpoonup}\bm{\nu}$* as $j\to\infty$, that is, \bigl{\langle}\nu^{k}_{n_{j}},g\bigr{\rangle}_{\mathcal{H}^{k}}\to\bigl{\langle}\nu^{k},g\bigr{\rangle}_{\mathcal{H}^{k}} for every $g\in\mathcal{H}^{k}_{0}([0,T),D;U)$ and every $k\in\mathbb{N}$* 2. (ii)

\bigl{\langle}\nu^{1}_{t},|\xi|^{p}\bigr{\rangle}_{\mathcal{H}^{1}}\leqslant c^{p}* for a.e. $t\in[0,T)$* 3. (iii)

\int_{0}^{T^{\prime}}\omega_{r}^{p}\big{(}\nu^{2}_{t}\big{)}\,dt\leqslant\liminf_{n\to\infty}\int_{0}^{T^{\prime}}\omega_{r}^{p}\big{(}\nu^{2}_{n,t}\big{)}\,dt* for every $r>0$ and $T^{\prime}\in[0,T)$* 4. (iv)

for $k\in\mathbb{N}$, let $\varphi\in L^{1}_{\textrm{loc}}([0,T)\times D^{k})$ and $\kappa\in C(U^{k})$ be nonnegative, and let $g(t,x,\xi):=\varphi(t,x)\kappa(\xi)$. Then

[TABLE] 5. (v)

Assume moreover that $D\subset\mathbb{R}^{d}$ is bounded, $T<\infty$ and that $\bm{\nu}_{n}$ have uniformly bounded support, in the sense that

[TABLE]

with $\mu_{t}^{n}\in\mathcal{P}_{T}(L^{p}(D;U))$ being dual to $\bm{\nu}_{n}$, then the following observables converge strongly:

[TABLE]

for every $g\in\mathcal{H}^{k,p}_{1}([0,T),D;U)$.

We skip the proof of this theorem as is very similar to that of Theorem 2.13.

Remark 2.22**.**

A closer look at the convergence statement (2.27) reveals that we can expect pointwise a.e convergence in space of the ensemble averages of the observable $g\in\mathcal{H}^{k,p}_{1}([0,T),D;U)$. On the other hand, time averaging in (2.27) seems essential. In other words, we have convergence of time averages of ensemble averages of the observables.

2.4. Statistical solutions

Using correlation measures we can now define statistical solutions of (1.1a). We need the following assumptions on the flux function in (1.1a),

[TABLE]

for some constant $C>0$ and $1\leqslant p<\infty$. The value of $p$ is given by available a priori bounds for solutions of (1.1a), for instance, from the entropy condition (1.3) (cf. (1.4)). For example, both the shallow water equations and the isentropic Euler equations are $L^{2}$-bounded, at least for solutions away from vacuum [10].

Statistical solutions are correlation measures (or equivalently, probability measures over $L^{p}$) satisfying the differential equation (1.1a) in a certain averaged sense. The full derivation can be found in [16], and we only provide the definition here.

Definition 2.23**.**

Let \bar{\mu}\in\mathcal{P}\big{(}L^{p}\big{(}D;U)\big{)} have bounded support,

[TABLE]

for some $M>0$. A statistical solution of (1.1a) with initial data $\bar{\mu}$ is a time-dependent map $\mu:[0,T)\mapsto\mathcal{P}(L^{p}(D;U))$ such that each $\mu_{t}$ has bounded support, and such that the corresponding correlation measures $(\nu^{k}_{t})_{k\in\mathbb{N}}$ satisfy

[TABLE]

in the sense of distributions, i.e.,

[TABLE]

for every \varphi\in C_{c}^{\infty}\big{(}D^{k}\times\mathbb{R}_{+},\ U^{\otimes k}\big{)} and for every $k\in\mathbb{N}$. (Here, $\bar{\bm{\nu}}$ denotes the correlation measure associated with the initial probability measure $\bar{\mu}$.)

Remark 2.24**.**

If the initial data $\bar{\mu}$ and a resulting statistical solution $\mu_{t}$ are both atomic, i.e. $\bar{\mu}=\delta_{\bar{u}}$ and $\mu_{t}=\delta_{u}$ with $\bar{u}\in L^{p}(D;U)$ and $u\in L^{p}((0,T)\times D;U)$, then it is easy to see that a statistical solution in the above sense reduces to a weak solution of (1.1a). Thus, weak solutions are statistical solutions.

Remark 2.25**.**

The evolution equation for the first correlation marginal of the statistical solution, i.e., for $k=1$ in (2.30), is equivalent to the definition of a measure-valued solution of (1.1a) [12, 15]. Thus, a statistical solution can be thought of as a measure-valued solution augmented with information about all possible multi-point correlations. Hence, a priori, a statistical solution contains significantly more information than a measure-valued solution.

3. Dissipative Statistical solutions and weak-strong uniqueness

In analogy with weak solutions, it is necessary to impose additional admissibility criteria for statistical solutions in order to ensure uniqueness and stability. In [16], the authors proposed an entropy condition for statistical solutions of scalar conservation laws. This condition was based on a non-trivial generalization of the Kruzkhov entropy condition to the framework of time-parameterized probability measures on $L^{1}(D)$. It was shown in [16] that these entropy statistical solutions were unique and stable in the 1-Wasserstein metric on $\mathcal{P}(L^{1}(D))$, with respect to perturbations of the initial data.

Although one can extend the entropy condition of [16] to statistical solutions for systems of conservation laws (1.1a), it is not possible to obtain uniqueness and stability of such entropy statistical solutions. Instead, one has to seek alternative notions of stability for systems of conservation laws.

A possible weaker framework for uniqueness (stability) is that of weak-strong uniqueness, see [47, 10] and references therein. Within this framework, one imposes certain entropy conditions and proves that the resulting entropy solutions will coincide with a strong (classical) solution if such a solution exists. Weak-strong uniqueness for systems of conservation laws with strictly convex entropy functions is shown in [10]. In fact, one can even prove weak-strong uniqueness results for the much weaker notion of entropy or dissipative measure-valued solutions of systems of conservation laws, see [41, 6, 15].

Our aim in this section is to propose a suitable notion of dissipative statistical solutions and prove a weak-strong uniqueness result for such solutions. Stability of solutions will be measured in the Wasserstein distance, whose definition we recall first.

Definition 3.1**.**

Let $X$ be a separable Banach space and let $\mu,\rho\in\mathcal{P}(X)$ have finite $p$th moments, i.e. $\int_{X}|x|^{p}d\mu(x)<\infty$ and $\int_{X}|x|^{p}d\rho(x)<\infty$. The $p$-Wasserstein distance between $\mu$ and $\rho$ is defined as

[TABLE]

where the infimum is taken over the set $\Pi(\mu,\rho)\subset\mathcal{P}(X^{2})$ of all transport plans from $\mu$ to $\rho$, i.e. those $\pi\in\mathcal{P}(X^{2})$ satisfying

[TABLE]

(see e.g. [45]).

As in [16, Section 4], our entropy condition for statistical solutions will rely on a comparison with probability measures that are convex combinations of Dirac masses, i.e. $\rho\in\mathcal{P}(L^{2}(D))$ such that $\rho=\sum_{i=1}^{M}\alpha_{i}\delta_{u_{i}}$ for coefficients $\alpha_{i}\geqslant 0$, $\sum_{i}\alpha_{i}=1$ and functions $u_{1},\dots,u_{M}\in L^{2}(D)$. From [16, Lemma 4.2], we observe that whenever $\rho$ is of this $M$-atomic form, there is a one-to-one correspondence between transport plans $\pi\in\Pi(\mu,\rho)$ and elements of the set

[TABLE]

defined for any $\alpha=(\alpha_{1},\dots,\alpha_{M})\in\mathbb{R}^{M}$ satisfying $\alpha_{i}\geqslant 0$ and $\sum_{i=1}^{M}\alpha_{i}=1.$ The set $\Lambda(\alpha,\mu)$ is never empty since $(\mu,\dots,\mu)\in\Lambda(\alpha,\mu)$ for any choice of coefficients $\alpha_{1},\dots,\alpha_{M}$. Note that the set $\Lambda(\alpha,\mu)$ depends on the target measure $\rho$ only through the weights $\alpha_{1},\dots,\alpha_{M}$.

Using this decomposition of transport plans with respect to M-atomic probability measures, we define the notion of dissipative statistical solution as follows.

Definition 3.2**.**

Assume that the system of conservation laws (1.1a) is equipped with an entropy function $\eta$. A statistical solution $\mu$ of (1.1a) is a dissipative statistical solution if

1. (i)

for every choice of coefficients $\alpha_{1},\dots,\alpha_{M}>0$ with $\sum_{i=1}^{M}\alpha_{i}=1$ and for every $(\bar{\mu}_{1},\dots,\bar{\mu}_{M})\in\Lambda(\alpha,\bar{\mu})$, there exists a function $t\mapsto(\mu_{1,t},\dots,\mu_{M,t})\in\Lambda(\alpha,\mu_{t})$, such that each measure $\mu_{i}\in\mathcal{P}_{T}(L^{p}(D;U))$ is a statistical solution of (1.1a) with initial data $\bar{\mu}_{i}$, 2. (ii)

for all test functions $0\leqslant\theta(t)\in C_{c}^{\infty}(\mathbb{R}_{+})$,

[TABLE]

We remark that the first condition in the above definition demands that the decomposition of a statistical solution into the components $\mu_{i}$ is still consistent with the underlying conservation law (1.1a). On the other hand, the second condition (3.2) amounts to requiring that the total entropy of $\mu$ decreases in time.

First, we investigate the stability of a dissipative statistical solution of (1.1a) with respect to statistical solutions built from finitely many classical solutions of (1.1).

Lemma 3.3**.**

Let $T>0$, set $p=2$, assume that

[TABLE]

(where we denoted by $f^{\prime\prime}$ the Hessian of $f$, i.e. $(f^{\prime\prime}(u))_{ijk}=\partial_{u^{j}}\partial_{u^{k}}f^{i}(u)$, $i,j,k=1,\dots,N$), and assume that the conservation law (1.1a) is equipped with an entropy pair $(\eta,q)$ for which

[TABLE]

(where $\eta^{\prime\prime}$ denotes the Hessian matrix of $\eta$)* for $c,C>0$. Let $\mu\in\mathcal{P}_{T}(L^{2}(D;U))$ be a dissipative statistical solution of (1.1a), and for $t\in[0,T)$ let $\rho_{t}=\sum_{i=1}^{M}\alpha_{i}\delta_{v_{i}(t)}$ for coefficients $\alpha_{i}>0$, $\sum_{i=1}^{M}\alpha_{i}=1$, and classical solutions $v_{1},\dots,v_{M}\in W^{1,\infty}(D\times[0,T);U)$ of (1.1a). Then*

[TABLE]

where $C=C(R)\geqslant 0$ is a constant only depending on $R:=\max_{i=1,\dots,M}\|v_{i}\|_{W^{1,\infty}(D\times\mathbb{R}_{+},U)}.$

Proof.

It is straightforward to verify that $\rho_{t}$ as defined above is a statistical solution of (1.1a) with initial data $\bar{\rho}:=\sum_{i=1}^{M}\alpha_{i}\delta_{\bar{v}_{i}}$, where $\bar{v}_{i}:=v_{i}(0)$.

Let $\bar{\mu}^{\ast}=(\bar{\mu}_{1}^{\ast},\dots,\bar{\mu}_{M}^{\ast})\in\Lambda(\alpha,\bar{\mu})$ define a transport plan that minimizes the transport cost between $\bar{\mu}:=\mu_{0}$ and $\bar{\rho}$, that is,

[TABLE]

(Here and in the remainder of this proof, we denote $L^{2}=L^{2}(D;U)$.) As $\mu_{t}$ is a dissipative statistical solution, there exists a map t\mapsto\big{(}\mu^{\ast}_{1,t},\dots,\mu^{\ast}_{M,t}\big{)}\in\Lambda(\alpha,\mu_{t}) such that

[TABLE]

for every $\varphi_{1},\dots,\varphi_{M}\in C_{c}^{\infty}(D\times[0,T))$. For each $1\leqslant i\leqslant M$, we have that

[TABLE]

Multiplying the above with $\alpha_{i}$ and summing over $i$, we obtain

[TABLE]

Subtracting (3.8) from (3.7) and choosing as a test function the vector-valued function $\varphi_{i}=\eta^{\prime}(v_{i}(x,t))\theta(t)$ for some scalar test function $0\leqslant\theta(t)\in C_{c}^{\infty}(\mathbb{R}_{+})$ (here, $\eta^{\prime}$ denotes the vector-valued derivative of $\eta$ with respect to $u$), and using the fact that

[TABLE]

yields

[TABLE]

As $v_{i}$ is a classical solution of (1.1a) and ${\mu}_{1,t}^{\ast},\dots,{\mu}_{M,t}^{\ast}$ are probability measures, we obtain from the entropy conservation of $v_{i}$ that

[TABLE]

for the same test function $\theta$ used in (3.9).

Since $\mu_{t}$ is a dissipative statistical solution of (1.1a), we have

[TABLE]

Subtracting (3.10) and (3.9) from the above inequality we obtain

[TABLE]

As $\eta$ is strictly convex (the lower bound in (3.4)), we have

[TABLE]

Similarly, using (3.3) and the fact that $\|v_{i}\|_{W^{1,\infty}}\leqslant R$, we obtain that

[TABLE]

Using the above estimates in (3.11) and choosing the test function $\theta(s)=\chi_{(0,t]}(s)$ (by approximating with smooth functions) yields

[TABLE]

Applying the integral form of Grönwall’s inequality to the above estimate results in

[TABLE]

the last step following from (3.6). As $(\mu_{1,t}^{*},\dots,\mu_{M,t}^{*})\in\Lambda(\alpha,\mu_{t})$, the above inequality implies (3.5) and concludes the proof. ∎

The estimate (3.5) implies stability of dissipative statistical solutions with respect to probability measures that are convex combinations of Dirac masses, concentrated on classical solutions of (1.1a). We can extend such a stability result to a more general class of strong solutions:

Definition 3.4**.**

A statistical solution $\mu$ is a strong statistical solution of (1.1a) if there is some $R>0$ such that for every $n\in\mathbb{N}$, there exists a $\rho_{n}\in\mathcal{P}_{T}(L^{2}(D,U))$ of the form $\rho_{n,t}=\sum_{i=1}^{N_{n}}\alpha_{i}\delta_{v_{i}(t)}$ such that each $v_{i}$ is a classical solution of (1.1a) and $v_{i}\in B_{R}\subset W^{1,\infty}(D\times(0,T))$, for all $1\leqslant i\leqslant N_{n}$ and such that

[TABLE]

We can now prove our main weak-strong uniqueness result:

Theorem 3.5**.**

Let $\bar{\mu}\in\mathcal{P}(L^{2}(D;U))$. Then under the same assumptions as in Lemma 3.3, if there exists a strong statistical solution $\mu$ of (1.1a), then it is unique in the class of dissipative statistical solutions.

Proof.

For a fixed $n\in\mathbb{N}$, Definition 3.4 implies the existence of a statistical solution $\rho_{t}=\sum_{i=1}^{N_{n}}\alpha_{i}\delta_{v_{i}(t)}$ such that

[TABLE]

Moreover, $v_{i}(x,t)\in B_{R}$ is a classical solution of (1.1a) with initial data $\bar{v}_{i}$ for all $i=1,\dots,N_{n}$.

Let $\gamma_{t}$ be another dissipative statistical solution of (1.1a) with $\gamma_{0}=\bar{\mu}$. By Lemma 3.3 we have

[TABLE]

for some $C=C(R)$. Using (3.14) and the triangle inequality yields

[TABLE]

Letting $n\rightarrow\infty$ concludes the proof of uniqueness of strong solutions. ∎

Remark 3.6**.**

If we assume that the initial data $\bar{\mu}$ is such that ${\rm supp}(\bar{\mu})\subset B_{R_{0}}$, with $B_{R_{0}}$ being the ball of radius $R_{0}>0$ in $W^{1,\infty}(D;U)$, then by classical results on local well-posedness for (1.1a) with strictly convex entropies [10], there exist $T(R_{0}),R(R_{0})>0$, such that for every initial data $\bar{v}\in{\rm supp}(\bar{\mu})$, there exists a corresponding classical solution $v\in W^{1,\infty}(D\times[0,T(R_{0})])$ and $\|v\|_{W^{1,\infty}(D\times[0,T(R_{0})])}\leqslant R(R_{0})$. Moreover, the data-to-solution map $\mathcal{S}_{t}:\operatorname{supp}\bar{\mu}\to L^{2}(D;U)$, is well defined for all $0\leqslant t\leqslant T(R_{0})$, and continuous because

[TABLE]

*Letting $\mu_{t}=\mathcal{S}_{t}\#\bar{\mu}$ for all $0\leqslant t\leqslant T(R_{0})$, one can verify that $\mu_{t}$ is indeed a dissipative statistical solution of (1.1a). Moreover, it is strong in the sense of Definition 3.4. Consequently, we can establish that *as long as the underlying initial data is supported on smooth functions, the resulting statistical solutions are locally well-posed. **

4. Numerical approximation of statistical solutions.

In this section, we will construct statistical solutions for the system of conservation laws (1.1a) by proposing an algorithm to numerically approximate it. We show, under reasonable hypotheses on the underlying numerical schemes, that the approximations constructed by this algorithm converge to a (dissipative) statistical solution of (1.1a). As in [15, 19, 17], the algorithm will be based on a finite volume spatio-temporal discretization and Monte Carlo sampling of the underlying probability space. The spatial domain $D$ will everywhere be assumed to be bounded.

4.1. Multidimensional finite volume framework

In this section we briefly describe numerically approximating conservation laws with finite volume and finite difference methods. For a complete review, one can consult [26, 28, 29, 34].

We discretize the computational spatial domain as a collection of cells

[TABLE]

with corresponding cell midpoints

[TABLE]

For simplicity we assume that our mesh is equidistant, that is,

[TABLE]

for some ${\Delta}>0$. For each cell, marked by ${\mathbf{i}}=(i^{1},\ldots,i^{d})$, we let $u^{{\Delta}}_{\mathbf{i}}(t)$ (and equivalently $u^{{\Delta}}_{i^{1},\ldots,i^{d}}(t)$) denote the averaged value in the cell at time $t\geqslant 0$. We consider the following semi-discrete scheme,

[TABLE]

where ${\mathbf{e}}_{1},\dots,{\mathbf{e}}_{d}$ are the canonical unit vectors in $\mathbb{R}^{d}$, and $F^{m,{\Delta}}$ is a numerical flux function in direction $m=1,\dots,d$. We say that the scheme is a $(2q+1)$-point scheme, when the numerical flux function $F^{m,{\Delta}}$ can be written as a function of $u^{{\Delta}}_{{\mathbf{i}}+j{\mathbf{e}}_{m}}(t)$ for $j=-q+1,\ldots,q$. We furthermore assume the numerical flux function is consistent with $f$ and locally Lipschitz continuous, which amounts to requiring that for every bounded set $K\subset\mathbb{R}^{N}$, there exists a constant $C>0$ such that for $m=1,\ldots,d$,

[TABLE]

whenever $u^{{\Delta}}_{{\mathbf{i}}-(q-1){\mathbf{e}}_{m}}(t),\ldots,u^{{\Delta}}_{{\mathbf{i}}+q{\mathbf{e}}_{m}}(t)\in K$. For the sake of notational simplicity we will write

[TABLE]

We let $\mathcal{S}^{{\Delta}}_{t}:L^{p}(D)\to L^{p}(D)$ be the spatially discrete numerical evolution operator defined by (4.1), mapping $\bar{u}\mapsto u^{\Delta}(t)$. Since $\mathcal{S}^{{\Delta}}_{t}$ is the composition of a projection from $L^{p}$ onto piecewise constant functions and a continuous evolution under an ordinary differential equation, we see that $\mathcal{S}^{{\Delta}}_{t}$ is measurable.

The current form of (4.1) is continuous in time, and one needs to employ a time stepping method to discretize the ODE system in time, usually through some strong stability preserving Runge–Kutta method [27].

As the operator $\mathcal{S}^{{\Delta}}$ is a measurable map, we can define an approximation of a statistical solution of (1.1a) with initial data $\bar{\mu}$ by

[TABLE]

Henceforth, $\mu_{t}^{{\Delta}}$ is referred to as an approximate statistical solution.

4.2. Convergence of approximate statistical solutions.

In this section, we investigate the convergence of the approximate statistical solutions $\mu_{t}^{{\Delta}}$ as the mesh is refined, i.e., as ${\Delta}\rightarrow 0$.

Theorem 4.1**.**

Consider the system of conservation laws (1.1a) with initial data $\bar{\mu}\in\mathcal{P}(L^{p}(D;U))$ for some $1\leqslant p<\infty$, such that $\operatorname{supp}(\bar{\mu})\subset B_{R}(0)\subset L^{p}(D;U)$, with $B_{R}(0)$ being the ball of radius $R$ and center 0, for some $R>0$. Assume that the semi-discrete finite volume scheme (4.1) satisfies the following conditions:

• (i)

$L^{p}$ bounds*:*

[TABLE]

• (ii)

Weak BV bounds*: There exists $s\geqslant p$ such that*

[TABLE]

with the constant $C=C(\|\bar{u}\|_{p})$ only depending on the $L^{p}$-norm of the initial data $\bar{u}$.

• (iii)

Approximate scaling*: There exists a constant $C>0$, possibly depending on the initial data $\bar{\mu}$ but independent of the grid size ${\Delta}$, such that for every $\ell>1$*

[TABLE]

*Here, $S^{p}_{r}(\mu)$ is the structure function associated with the time parameterized probability measure $\mu_{t}\in\mathcal{P}_{T}(L^{p}(D;U))$ *(equivalently, time-dependent correlation measure $\bm{\nu}\in\mathcal{L}^{p}([0,T),D;U)$), defined as

[TABLE]

Then there is a subsequence ${\Delta}^{\prime}\to 0$ such that the approximate statistical solutions $\mu^{{\Delta}^{\prime}}$ converge strongly to some $\mu\in\mathcal{P}_{T}(L^{p}(D,U))$, in the sense of Theorem 2.21(v).

Proof.

We will show that that the approximate statistical solutions $\mu_{t}^{{\Delta}}$, defined in (4.3), satisfy the conditions of Theorem 2.21 and hence converge (up to a subsequence). To this end, we can readily verify from the uniform $L^{p}$ bounds (4.4) and the fact that $\operatorname{supp}\bar{\mu}\subset B_{R}(0)$, that

[TABLE]

with $C$ being the constant in (4.4).

Let $u^{\Delta}(t)=\mathcal{S}^{{\Delta}}_{t}\bar{u}$ be the solution generated by the scheme (4.1). Denoting

[TABLE]

we obtain from the weak BV estimate (4.5) and Hölder’s inequality that (recall that $D$ was assumed to be bounded)

[TABLE]

The above inequality holds for every $\bar{u}\in B_{R}(0)\subset L^{p}(D;U)$.

Next, for any $r\leqslant{\Delta}$, a straightforward but tedious calculation yields

[TABLE]

where $C_{d}=3^{d-1}$ results from successive applications of the triangle inequality. Hence, summarizing the above calculation, we obtain that for any $r\leqslant{\Delta}$,

[TABLE]

for a constant $C$ that depends on the dimension, the support of the initial probability measure and the final time but is independent of the grid size ${\Delta}$.

Now, for any $\ell>1$ and $r=\ell{\Delta}$, we have

[TABLE]

Here, the constant is independent of ${\Delta}$. By combining the above estimate with (4.10), we obtain that

[TABLE]

for any $r>0$.

Given the independence of the constant in (4.11) with respect to the grid size ${\Delta}$, we see from (4.11) that the condition of uniform diagonal continuity (2.24) is satisfied. Hence, up to a subsequence still indexed by ${\Delta}$, $\mu^{{\Delta}}_{t}$ converges to some $\mu\in\mathcal{P}_{T}(L^{p}(D;U))$. ∎

Several remarks on the assumptions in the above theorem follow.

Remark 4.2**.**

There are many examples of finite volume/difference schemes of the form (4.1) which satisfy the uniform $L^{p}$ bound (4.4) and the weak BV bound (4.5). Assume that the system of conservation laws (1.1a) possesses an entropy function $\eta$ that satisfies

[TABLE]

for constants $C_{1},C_{2}>0$ and $p\in[1,\infty)$. Then the uniform $L^{p}$ bound (4.4) follows for any scheme of form (4.1) that satisfies a discrete entropy inequality,

[TABLE]

with a numerical entropy flux $Q^{m,{\Delta}}$ that is consistent with the entropy flux $q^{m}$ in (1.3) for $1\leqslant m\leqslant d$.

In many cases the weak-BV bound (4.5) also follows from the discrete entropy inequality (4.13). Examples of schemes which satisfy the discrete entropy inequality (4.13) and the weak BV bound (4.5) are the so-called entropy stable Lax–Wendroff schemes and the TeCNO schemes of [18].

Remark 4.3**.**

The approximate scaling assumption (4.6) can be thought of as a weaker version of the so-called self-similarity at small scales assumption of Kolmogorov in his K41 theory for fully developed turbulence in incompressible fluid flows, see hypothesis H2, equation (6.3), page 75 of [24]. Kolmogorov based his hypothesis on the fact that solutions of the incompressible Navier–Stokes (Euler) equations scale exactly. Similar considerations also apply to several prototypical examples of systems of conservation laws (1.1a). In particular, for the compressible Euler equations (5.1) (in any space dimension), it can be readily checked that if $u(x,t)$ is a solution, then $\ell^{\theta}u(\ell x,\ell t)$ is also a solution for any $\theta,\ell>0$. Hence, it is reasonable to hypothesize scaling, analogous to the Kolmogorov hypothesis, for systems of conservation laws.

It is essential to also point out the differences in our hypothesis (4.6) to the standard Kolmogorov hypothesis for turbulent incompressible flows. First, our hypothesis pertains only to the numerical solution, generated by the finite volume scheme (4.1). Moreover, we require mere inequalities in the scaling law (4.6), in contrast to the standard Kolmogorov hypothesis of equality.

Remark 4.4**.**

Intermittency is widely accepted to be a characteristic of turbulent flows, see [24]. It is believed that intermittency stems from the fact that turbulent solutions do not scale exactly as in the Kolmogorov hypothesis. We automatically incorporate a form of intermittency by only requiring an upper bound in (4.6), instead of an equality. Hence, the scaling exponent in (4.6) can depend explicitly on the underlying length scale, provided that it is bounded above by $1/s$. This encodes a form of intermittency in the approximate solutions.

Remark 4.5**.**

Another approach to incorporating intermittency and relaxing the scaling condition (4.6) is to consider a decomposition of the approximate statistical solution $\mu^{{\Delta}}_{t}$ into a mean flow and a fluctuation. Defining the mean flow by

[TABLE]

we see that the mean flow is well defined for almost every $(x,t)\in D\times(0,T)$. Similarly, we can define fluctuations of $\tilde{\mu}^{\Delta}\in\mathcal{P}_{T}(L^{p}(D;U))$ by its action on all observables $g\in\mathcal{H}^{k,p}([0,T],D;U)$,

[TABLE]

We can relax the assumption (4.6) by requiring that only the structure function associated with the fluctuation scales approximately, i.e.,

[TABLE]

If we further assume that the mean flow is BV and $L^{\infty}$, i.e.,

[TABLE]

for some constant that is independent of the mesh size ${\Delta}$, then a straightforward but tedious calculation yields for any $r=\ell\Delta$

[TABLE]

Thus, the condition (2.24) in Theorem 2.21 is satisfied in this case. A similar argument can be made by imposing some form of (uniform) Hölder continuity on the mean flow.

4.3. Consistency of the numerical method

We fix an initial measure $\bar{\mu}\in\mathcal{P}(L^{p}(D,U))$. For any $u\in L^{p}(D,U)$ we define the local average of $u$ as

[TABLE]

where $|\mathcal{C}_{\mathbf{i}}|$ denotes the $d$-dimensional Lebesgue measure of $\mathcal{C}_{\mathbf{i}}$. We now state the “Lax–Wendroff theorem” for our numerical method, that is, consistency of the method with the PDE.

Theorem 4.6** (Lax–Wendroff theorem for statistical solutions).**

Let the initial data $\bar{\mu}$ have bounded support, $\operatorname{supp}\bar{\mu}\subset B_{K}(0)\subset L^{p}(D,U)$ for some $K>0$. Let $\mu^{\Delta}$ be given by (4.3) for ${\Delta}>0$, and assume that for some sequence $\Delta_{n}\to 0$, the sequence $\{\mu^{\Delta_{n}}\}_{n\in\mathbb{N}}$ converges strongly to $\mu$ in $\mathcal{P}_{T}(L^{p}(D,U))$, in the sense of Theorem 2.21(v). Assume moreover that the following weak BV bound is fulfilled:

[TABLE]

for some $0<s$. Then $\mu_{t}$ is a statistical solution of (1.1).

Given the complicated notation and very technical nature of the proof of the above theorem, we illustrate the main steps of the proof in a very special case, namely $k=2$ for a one-dimensional scalar conservation law ($d=N=1$). The proof in the general case is postponed to Appendix D.

Proof for the second moment of a scalar conservation law in one spatial dimension.

We consider a scalar conservation law ($N=1$) in one spatial dimension ($d=1$). By (4.4), there is some $K_{T}>0$ such that

[TABLE]

Let $\{u_{i}^{\Delta}\}_{i\in\mathbb{Z}}$ be computed by (4.1). Denote $F^{\Delta}=F^{1,{\Delta}}$ and, for $u\in L^{p}(D)$, write

[TABLE]

For all pairs of cells $i,j$, we have, by the product rule,

[TABLE]

Hence, for arbitrary $\varphi\in C_{c}^{\infty}(\mathbb{R}^{2}\times[0,T))$, we get

[TABLE]

Multiply by ${\Delta}^{2}$, sum over all $i,j\in\mathbb{Z}$ and perform a summation-by-parts to obtain

[TABLE]

where we have denoted $\partial_{1}^{\Delta}\varphi(x,y,t)=\frac{\varphi(x+{\Delta},y,t)-\varphi(x,y,t)}{{\Delta}}$, and similarly for $\partial_{2}^{\Delta}\varphi$. From the special form (4.3) of $\mu_{t}^{\Delta}$, we therefore have

[TABLE]

We write now

[TABLE]

The last term vanishes as ${\Delta}\to 0$, since

[TABLE]

A similar computation holds for the integral involving $F_{j+{\nicefrac{{1}}{{2}}}}^{\Delta}(u)$. Setting ${\Delta}={\Delta}_{n}$ then gives

[TABLE]

which completes the proof. ∎

Remark 4.7**.**

We can readily show that the limit statistical solution $\mu_{t}$ is a dissipative statistical solution, assuming that the underlying finite volume method satisfies the discrete entropy inequality (4.13). To this end, for every choice of coefficients $\alpha_{1},\dots,\alpha_{M}>0$ with $\sum_{i}^{M}\alpha_{i}=1$ and every $\left(\bar{\mu}_{1},\ldots,\bar{\mu}_{M}\right)\in\Lambda(\alpha,\bar{\mu})$, we construct $\mu^{{\Delta}}_{i,t}=\mathcal{S}^{{\Delta}}_{t}\#\bar{\mu}_{i}$ as the approximate statistical solution generated by the scheme (4.1). By the convergence theorem 4.1, we can show that each $\mu^{{\Delta}}_{i,t}$ converges (possibly along a further subsequence), in the topology of Theorem 2.21, to $\mu_{i,t}\in\mathcal{P}_{T}(L^{p}(D;U))$ as ${\Delta}\rightarrow 0$. By Theorem 4.6, each $\mu_{i}$ is a statistical solution of (1.1a) with initial data $\bar{\mu}_{i}$, and the condition (3.2) is a straightforward consequence of the discrete entropy inequality (4.13) and the growth condition (4.12) with $p=2$.

4.4. Monte Carlo algorithm

While (4.3) provides an abstract definition of $\mu_{t}^{{\Delta}}$, it is not amenable to practical computations, since it requires the computation of the trajectory of the numerical solution operator for almost all possible initial data $\bar{u}\in\operatorname{supp}\bar{\mu}$. We will further approximate $\mu_{t}^{{\Delta}}$ by sampling it for a large ensemble of initial data, drawn from the initial probability measure.

The Monte Carlo algorithm has been shown to be robust in tackling high dimensional problems with low regularity [39, 40], and has later been demonstrated to perform very well for computing measure valued solutions [15].

Algorithm 1** (Monte Carlo Algorithm).**

In the rest of this paper we will refer to the above algorithm simply as “the Monte Carlo Algorithm”. Note that for any admissible observable $g\in\mathcal{H}^{k,p}_{1}([0,T],D;U)$, using (4.3) and (4.20), we obtain that

[TABLE]

Remark 4.8**.**

One should note that the probability measure $\mu^{\Delta,M}_{T}$ is indeed a random probability measure depending on some probability space $\Omega$ from which $\bar{u}_{1},\ldots\bar{u}_{M}$ are being drawn.

Using well known results for weak convergence of Monte Carlo [44], we can prove that the Monte Carlo approximation of the statistical solution converges as the number of samples is increased.

Theorem 4.9**.**

Let $\bar{\mu}\in\mathcal{P}(L^{p}(D,U))$ have bounded support, let $\mathcal{S}^{{\Delta}}$ be some numerical evolution operator, and let $\mu_{t}^{{\Delta},M}$ be defined through the Monte Carlo Algorithm. Let $\mu_{t}^{{\Delta}}$ be defined by (4.3). Then for every admissible observable $g\in\mathcal{H}^{k,p}_{1}([0,T],D;U)$, we have

[TABLE]

where

[TABLE]

The proof of the above theorem follows by standard arguments for proving convergence of Monte Carlo approximations. It is analogous to the proof of Monte Carlo convergence to statistical solutions of scalar conservation laws, see Theorem 2 of [17]. Note that the right hand side in (4.22) is bounded on account of the hypothesis (2.16) on admissible observables $g\in\mathcal{H}^{k,p}_{1}([0,T],D;U)$.

5. Numerical experiments

For all the numerical experiments in this section, we consider the two-dimensional compressible Euler equations,

[TABLE]

The system is closed with the equation of state

[TABLE]

We set $\gamma=1.4$ for all experiments.

5.1. Kelvin–Helmholtz problem

We start with this well-known test case for the development of instabilities in fluid flows, which was also extensively studied in [15].

The Kelvin–Helmholtz initial data is a shear flow, separating two states of varying density and pressure,

[TABLE]

We assign periodic boundary conditions, and the two states are given as $\rho_{L}=2$, $\rho_{R}=1$, $w^{x}_{L}=-0.5$, $w^{x}_{R}=0.5$, $w^{y}_{L}=w^{x}_{L}=0$ and $p_{L}=p_{R}=2.5$. The interfaces between the two states are given as

[TABLE]

where $K=10$, $\epsilon>0$, and $\{a^{i}_{j}\}$ and $\{b^{i}_{j}\}$ are uniformly distributed random variables on the interval $[0,1]$. We normalize the $a^{i}_{j}$ such that $\sum_{j}a^{i}_{j}=1$.

The initial measure $\bar{\mu}$ is given by the distribution of the random field $u_{0}$. Note that although $\bar{\mu}$ is a probability measure on the infinite-dimensional space $L^{p}(D;U)$, it is only concentrated on a 40-dimensional subset of this space.

As was already shown in [15], there is no convergence for single realizations (samples) of the problem (5.1),(5.2) as the mesh is refined. We observe this behaviour from Figure 1, where we display the approximate density at time $T=2$, computed with a second-order high-resolution finite volume scheme using an HLLC approximate Riemann solver and WENO reconstruction, together with a second-order SSP Runge–Kutta time integrator, on a sequence of successively refined meshes. As seen from the figure, structures at finer and finer scales are generated upon mesh refinement, impeding convergence. This lack of convergence is also verified from Figure 2 (A), where the so-called Cauchy rates of the density, i.e., quantities of the form

[TABLE]

with $\Psi^{\Delta}$ being any function computed with (4.1) on a mesh with mesh size ${\Delta}$.

On the other hand, the theory developed in Section 4.1 suggests that observables $g\in\mathcal{H}^{k,p}_{1}([0,T],D;U)$, for all $k$, should converge on mesh refinement. We start by considering observables with respect of the first marginal $\nu^{1}$ of the underlying approximate statistical solution. In particular, we consider the mean and variance given by

[TABLE]

The above quantities are defined a.e. in $D\times[0,T]$ and $\nu^{1,{\Delta}}$ is the first marginal of the approximate statistical solution $\mu_{t}^{{\Delta},M}$ generated by the Monte Carlo Algorithm. It is straightforward to check that the mean and the variance are admissible observables, in the sense of Convergence Theorem 2.21.

We plot the mean and the variance of the density at time $T=2$ in Figure 3. As seen from this figure, and in contrast to single samples, the mean and variance clearly converge upon mesh refinement. Moreover, the variance is also concentrated along the so-called mixing zone, which spreads out from the initial interface. The convergence of the mean and the variance is further verified from Figure 2(B,C), where the Cauchy rates (5.4) of the mean and variance $M^{\Delta},V^{\Delta}$ are displayed as a function of mesh resolution.

To quantify the convergence of the distribution of $\mu^{{\Delta},M}$ we consider the following Cauchy rates,

[TABLE]

with $W_{p}$ being the Wasserstein metric defined in (3.1) and $\nu^{k,{\Delta}}$, the $k$-th correlation marginal, corresponding to the (approximate) statistical solution $\mu_{t}^{{\Delta},M}$ generated by the Monte Carlo Algorithm. One can check that for all $k\in\mathbb{N}$

[TABLE]

(see Appendix E for a proof). As the Wasserstein metric metrizes the weak topology on probability measures, we may conclude from Theorems 2.21 and 4.1 that under the assumptions of some form of time continuity, the right hand side of (5.7) goes to zero as ${\Delta}\rightarrow 0$. This convergence is verified in Figure 4 (A), where we plot the corresponding Cauchy rates for the distance (5.6) with respect to the density (see Appendix F for details about how the Wasserstein distance was computed numerically).

Next, we consider computation of observables with respect to the second correlation marginal $\nu^{2,\Delta}_{x,y,t}$ of the approximate statistical solution $\mu^{{\Delta},M}_{t}$. The most interesting observable in this regard is the approximate structure function $S^{p}_{r}(\mu^{{\Delta},M}_{t},T)$ (4.7). This is clearly an admissible observable in the sense of Theorem 2.21. For computational purposes, it is easier to compute the time-sections of the structure function, namely

[TABLE]

for $t\in[0,T)$. In Figure 5 we plot \omega^{2}_{r}\big{(}\nu^{2,{\Delta},M}_{T}\big{)}^{1/p} for a sequence of mesh sizes ${\Delta}$. Moreover, we consider three different setups in the figure. In Figure 5 (A) and (B), we set $\epsilon=0.1$ in (5.3) and $T=2$ and $T=4$, respectively, and in Figure 5 (C) we set $\epsilon=0.01$ and $T=2$. As seen from the figures, it is clear that the approximate structure functions converge as the mesh is refined. Moreover, the structure functions (approximately) behave as

[TABLE]

The computed values of $\theta_{p}(T)$ are seen in the legend in Figure 5.

The numerical convergence of structure functions is further verified in Figure 4 (B), where we plot the Cauchy rates (5.4) with g(r,T)=\omega^{2}_{r}\big{(}\nu^{2,{\Delta},M}_{T}\big{)} with $p=1$. In this figure $T=2$ and $\epsilon=0.01$. Clearly, the structure function (as a function of length scale $r$) converges as the mesh is refined.

Finally, we compute the Wasserstein Cauchy rates ${\mathcal{W}}_{1}^{2,\Delta}(T)$ (5.6) with respect to the density, over successively refined mesh sizes and display the result in figure 4 (c). This figure clearly shows that there is convergence (with respect to mesh resolution and number of Monte Carlo samples) in this metric.

We are interested in computing the two-point correlation marginal $\nu^{2,{\Delta},M}_{x,y,t}$ for point pairs $(x,y)\in D$. In this context, we realize this Young measure by empirical histograms, plotted in Figure 6. In this figure we show the empirical histogram of the two-point correlation Young measure of the density, on successively refined grids, for two different point pairs $x=(0.7,0.7),$ $y=(0.4,0.2)$ (top 2 rows) and $x=(0.7,0.7)$, $y=(0.7,0.8)$ (bottom 2 rows). We see from this figure that the empirical histograms of the two-point correlation marginals converge (visually) on mesh refinement. Moreover, there is a clear difference in the correlation structures at different point pairs.

5.2. Richtmeyer–Meshkov problem

Our second test case is the well-studied Richtmeyer–Meshkov problem (see [15] and references therein), which involves a very complicated solution of the compressible Euler equations (5.1), modeling the complex interaction of strong shocks with unstable interfaces. The underlying initial data is given as

[TABLE]

We assign periodic boundary conditions on $D=[0,1]^{2}$. The interface between the two states is given as

[TABLE]

where $K=10$, $\epsilon>0$, and $\{a_{j}\}$ and $\{b_{j}\}$ are uniform random variables on the interval $[0,1]$. We normalize the $a_{j}$ such that $\sum_{j}a_{j}=1$. The initial probability measure $\bar{\mu}$ is given by the law of the above random field, and lies in $\mathcal{P}(L^{p}(D))$ for every $1\leqslant p<\infty$.

As in the case of the Kelvin–Helmholtz problem, there is no convergence (on mesh refinement) for single samples (realizations). This non-convergence is demonstrated in Figure 7(A), where the Cauchy rates (5.4) with respect to the density at time $T=5$ are shown. We visualize the density for different mesh resolutions in Figure 8 (top row). As seen from this figure, the solution at this time is very complicated on account of the interaction between the incoming strong shock (which has been reflected, due to periodic boundary conditions) and the unstable interface, which leads to the generation of turbulent small scale eddies.

On the other hand, and as predicted by the convergence theory developed in Section 4.1, statistical observables such as the mean and the variance (5.5) converge on mesh refinement, as shown in Figure 8 (middle and bottom). Furthermore, this figure shows how the small scale structures are averaged out in the mean, whereas the small scale information is encoded in the variance, which is concentrated around the mixing zone. We also verify the convergence of the Wasserstein distance ${\mathcal{W}}_{1}^{2,\Delta}(T)$ with respect to the density for successively refined meshes in Figure 7 (B).

Next, we compute the time sections of the structure function $\omega_{r}^{2}(\nu^{2,{\Delta},M}_{t})$ defined in (5.8). These are shown in Figure 9, where we have used $T=5$, $r\in\left[1/1024,1/32\right]$, and two different values of the perturbation parameter in (5.10). As seen from the figure, the structure function clearly converges on mesh refinement. Moreover, it behaves as in (5.9), with exponents shown in Figure 9. The convergence of the structure function is further verified by plotting the Cauchy rates for the structure function, as a function of the length scale $r$ in Figure 7 (C).

In Figure 7(D) we plot the Wasserstein distance ${\mathcal{W}}_{1}^{2,\Delta}(T)$ for the density and $T=5$, on a sequence of successively refined meshes. As shown in the figure, this distance converges on mesh (and sample) refinement.

Finally, in Figure 10, we plot histograms that represent the two-point correlation measure for the density at two different point pairs and at time $T=5$. These histograms show that the two-point correlation structure for this initial datum is very different from the correlation structure for the Kelvin–Helmholtz problem (Figure 6).

5.3. Fractional Brownian motion

The initial probability measure $\bar{\mu}\in\mathcal{P}(L^{p}(D))$ in the previous two numerical experiments was realized as a probability measure on high- but finite-dimensional subsets. We now consider initial probability measures that are concentrated on genuinely infinite-dimensional subsets of $L^{p}(D)$.

We will assume that the initial probability measure for the two-dimensional compressible Euler equations (5.1) will correspond to a fractional Brownian motion. Introduced by Mandelbrot et al. [38], fractional Brownian motion can be seen as a generalization of standard Brownian motion with a scaling exponent different from $1/2$.

We consider the following initial data,

[TABLE]

where $B^{H}_{1}$ and $B^{H}_{2}$ are two independent two dimensional fractional Brownian motions with Hurst index $H\in(0,1)$. Standard Brownian motion corresponds to a Hurst index of $H=1/2$.

To generate fractional Brownian motion, we use the random midpoint displacement method originally introduced by Lévy [35] for Brownian motion, and later adapted for fractional Brownian motion [23, 46]. Consider the uniform grid $0=x_{\nicefrac{{1}}{{2}}}<\dots<x_{N+{\nicefrac{{1}}{{2}}}}=1$ with $x_{i+{\nicefrac{{1}}{{2}}}}=i{\Delta}$ and ${\Delta}=\frac{1}{N}$, where $N=2^{k}+1$ is the number of cells for some $k\in\mathbb{N}$. We first fix the corners

[TABLE]

Recursively update the values on the edges as

[TABLE]

and correspondingly for $w^{x,H,\Delta x}_{2^{k-l}j,2^{k-l-1}(2i+1)}(\omega;0)$. For the values in the center of the cells, we use the following expression:

[TABLE]

for $l=0,\ldots,k$ and for $i,j=0,\ldots,2^{l}$. Here $(X_{k,n})_{(k,n)\in\mathbb{N}^{2}}$ is a collection of normally distributed random variables with mean 0 and variance 1. That is, we bisect every cell and set the middle value to the average of the neighbouring values plus some Gaussian random variable. The same procedure is repeated for $w^{y,H,{\Delta}}$. See Figure 11 for a sample of the initial velocity field with standard Brownian motion, i.e., with $H=0.5$.

The initial probability measure is given by the law of the above random field and the dimension of its support increases with decreasing mesh size. Hence, in the limit ${\Delta}\to 0$ we are approximating a probability measure supported on an infinite dimensional subspace of $L^{2}(D;U)$

We compute the statistical solutions with Algorithm 1, with the fractional Brownian motion initial data for two different Hurst indices, $H=0.1$ and $H=0.5$. Statistical observables corresponding to the one-point correlation marginal, such as the mean and variance, converge on mesh refinement, as shown in Figure 12. We also plot the mean and variance of the density at the highest mesh resolution of $1024^{2}$ and time $T=0.25$, for the two different Hurst indices, in Figure 13. As seen from the figure, there is a clear difference in the spatial structure of the mean and the variance as the Hurst index is changed. Moreover, the spatial structure of these statistical quantities is much more complicated than in the case of the Kelvin–Helmholtz and Richtmeyer–Meshkov initial data, with no clear large scale structures such as shocks. On the other hand, the statistical quantities have more small-scale structures. This is more pronounced in the $H=0.1$ case than for standard Brownian motion.

For $r\in\left[1/1024,1/32\right]$, we plot the (time sections of) the structure function $\omega_{r}^{2}(\nu^{2,{\Delta},M}_{t})$ at $t=0.25$ in Figure 14. The structure functions clearly converge on mesh refinement. This is also verified for both Hurst indices in Figure 15 (A,B), where we plot the Cauchy rates (5.4) for the structure function, with respect to the length scale $r$. Moreover, the structure functions scale as in (5.9).

In Figure 15 (C,D), we plot the Cauchy rates with respect to the Wasserstein distance ${\mathcal{W}}^{2,{\Delta}}_{1}(0.25)$, with respect to grid resolution, for both Hurst indices. We verify from this figure that these distances also converge on mesh refinement and sample augmentation.

Finally, in Figure 16, we plot histograms representing the two-point correlation marginals of the density, computed on the finest grid resolution of $1024^{2}$, for two different point pairs. The figure shows that the two-point correlation structure is again very different for different Hurst indices, and from the correlation structures for the previous numerical experiments.

5.4. Stability of the computed statistical solution

A priori, the computed statistical solution depends on the specifics of the underlying initial probability measure $\bar{\mu}$ as well as on the details of the numerical scheme (4.1), used within the Monte Carlo Algorithm. We investigate the stability of the computed statistical solution with respect to these parameters in the specific case of the Kelvin–Helmholtz problem (5.2).

In [15], the authors had already demonstrated the stability of the computed measure-valued solution with respect to the variations of the underlying numerical method, or to the size and type of perturbations to the Kelvin–Helmholtz initial data. As the computed measure-valued solution in [15] is identical to the first correlation marginal of our computed statistical solution, we can assume that the observables with respect to the first correlation marginal are also stable. Therefore, we investigate the stability of observables with respect to the second correlation marginal $\nu^{2,{\Delta},M}$. The results are summarized below.

• •

Stability with respect to amplitude of perturbations. We vary the size of the perturbation parameter $\epsilon$ in (5.3) over two orders of magnitude, from $\epsilon=0.001$ to $\epsilon=0.1$. The computed (time section of) structure function (5.8) for $p=2$ and at time $T=2$, on the finest resolution of $1024^{2}$ points and $1024$ Monte Carlo samples is shown in Figure 17 (A, left). As seen from the figure, the computed structure functions are very close to each other and scale as (5.9) with $\theta\approx 0.61$. This indicates stability of the computed structure function with respect to the amplitude of perturbations in the initial data. This stability is further verified in Figure 17 (A, right) where we plot the Wasserstein distance \bigl{\|}W_{1}\left(\nu^{2,\epsilon}_{x,T},\nu^{2,\frac{\epsilon}{2}}_{x,T}\right)\bigr{\|}_{L^{1}(D^{2})} with respect to the density, at time $T=2$ for different values of the perturbation parameter. The plot shows (linear) convergence with the decay of the perturbation, indicating stability of the computed statistical solution vis a vis perturbation amplitude.

• •

Stability with respect to type of perturbations. In all the numerical experiments for the Kelvin–Helmholtz initial data, we have assumed that the random variables $a_{j},b_{j}$ in (5.3) are chosen from a uniform distribution. Here, we choose these random variables from a standard normal distribution. This amounts to varying the corresponding initial probability measure for (5.2). The consequent change in the structure function (5.8) for two different amplitudes of the perturbation parameter $\epsilon$ in (5.3) are shown in Figure 17 (B, left). The figure clearly shows that the computed structure functions are very close to the ones computed with the uniform distribution. This stability with respect to the type of perturbation is further verified in Figure 17 (B, right) where we plot the \big{\|}W_{1}\left(\nu^{2,\epsilon}_{x,T},\hat{\nu}^{2,\epsilon}_{x,T}\right)\big{\|}_{L^{1}(D^{2})} at time $T=2$. Here, $\nu,\hat{\nu}$ refer to the correlation measures, computed with the uniform and standard normal random variables, respectively. The plot shows convergence with the decay of the perturbation, indicating stability of the computed statistical solution, vis a vis perturbation type.

• •

Stability with respect to choice of numerical scheme. In order to investigate the stability of the computed statistical solutions to the choice of the underlying numerical scheme in (4.1), we vary the reconstruction procedure, i.e., we use a high-resolution finite volume scheme based on the HLLC flux, but with either MC or WENO2 reconstructions (see e.g. [34]). The choice of the reconstruction leads to change in the sub-grid scale numerical viscosity of the overall approximation. We plot the structure function (5.8) in Figure 17 (C, left) and observe a very minor change in the structure function. This issue is investigated further in Figure 17 (C, right) where the Wasserstein distances \big{\|}W_{1}\left(\nu^{2,{\Delta}}_{\text{WENO2},x,T},\nu^{2,{\Delta}}_{\text{MC},x,T}\right)\big{\|}_{L^{1}(D^{2})} at time $T=2$ are plotted. Here $\nu^{2,{\Delta}}_{\text{WENO2}}$ is the second correlation marginal, computed with the WENO2 reconstruction procedure, and $\nu^{2,{\Delta}}_{\text{MC}}$ is the second correlation marginal, computed with the MC reconstruction procedure. We observe convergence of this distance with respect to resolution. This allows us to conclude that the statistical solutions are stable with respect to the choice of the underlying numerical method, at least for this Kelvin–Helmholtz problem.

5.5. Statistical steady state and regularity

In addition to the four numerical experiments reported in the last section, i.e., Kelvin–Helmholtz, Richtmeyer–Meshkov, and fractional Brownian motion with two different Hurst indices of $H=0.1$ and $H=0.5$, we have performed two further numerical experiments. Both of them consider the two-dimensional compressible Euler equations with the following initial data:

• •

Fractional Brownian motion initial data with Hurst index $H=0.75$.

• •

Shock-vortex interaction initial data, see Section 6.3.2 of [18] and references therein.

The Monte Carlo Algorithm is used to compute the approximate statistical solution $\mu_{t}^{{\Delta},M}$ for these additional sets of initial data.

We focus on the (time-sections of) the structure function (5.8) and find that in all six numerical experiments, the structure function behaved as (5.9). The exponent $\theta_{p}(t)$ as a function of time, for $p=1,2,3$ and for each numerical experiment is shown in Figure 18. We observe the following from this figure.

• •

First, the exponent $\theta_{p}(t)$ reaches a steady state rather quickly, when compared to the dynamic behaviour of the solution. In other words, Figure 18 seems to suggest that statistical equilibrium is reached significantly faster than the (deterministic) steady state for individual realizations. Hence, the system evolves dynamically for each sample, while the whole ensemble has already reached statistical equilibrium. The time scale at which this statistical equilibrium is reached depends on the specifics of the initial data.

• •

For all the experiments except the shock-vortex interaction, there is a very interesting behavior of the structure function (5.8) with respect to time. In particular, the exponent $\theta_{p}(t)$ for $p=2,3$ clearly increases with time, indicating that the non-linear evolution statistically regularizes the solution in some manner. The exception is for the shock-vortex interaction where this exponent remains constant with time. This can be explained by the fact that the shock-vortex interaction results in a solution whose total variation (TV) norm is bounded. Hence, one can readily verify that $\theta_{p}(t)=\frac{1}{p}$, which is approximately realized in the computations. On the other hand and, as shown in Figure 19, the (average) BV-norm blows up for all the remaining test cases. Hence, in these problems, there is a complex interaction of structures at different length scales that leads to a subtle statistical regularity.

• •

Last, but not least, we observe that in all the numerical experiments we have considered, the structure functions scale as (5.9). Hence, the approximate scaling assumption (4.6) in Theorem 4.1 is always observed to be satisfied.

5.6. Reproducing the numerical experiments

All experiments were carried out using the open source Alsvinn simulator [36]. For a full description on how the experiments were carried out, along with the raw data and post processing scripts, consult

https://github.com/kjetil-lye/systemspaper_experiments

and

https://github.com/kjetil-lye/statistical_systems_paper_experiments

6. Discussion

We consider hyperbolic systems of conservation laws (1.1). Although the standard solution framework of entropy solutions has been shown to be well-posed for scalar conservation laws and one-dimensional systems, it is now clearly established that entropy solutions for multi-dimensional systems are not unique, nor are they amenable to numerical approximation. On the other hand, numerical evidence presented in [15, 19] and references therein suggests that a statistical notion of solutions might be more appropriate for (1.1), even if the initial data and other underlying parameters are deterministic.

Entropy measure-valued solutions [12, 15] are a possible solution framework. Although global existence and numerical approximation results for entropy measure-valued solutions are available, it is well-known that measure-valued solutions are not unique, even for scalar conservation laws. This is largely on account of lack of information about multi-point spatial correlations.

Inspired by the need to incorporate correlations, the authors of [16] proposed a framework of statistical solutions for hyperbolic systems of conservation laws. Statistical solutions are time-parameterized probability measures on the space of $p$-integrable functions. They were shown in [16] to be equivalent to adding information about all possible multi-point correlations to the measure-valued solution. The time-evolution of these measures is prescribed by a system of nonlinear tensorized moment-transport equations (2.30). Under an additional entropy condition, the well-posedness of statistical solutions for scalar conservation laws was shown in [16] and the numerical approximation of statistical solutions for scalar conservation laws was considered in [17].

Our main aim in this paper was to propose a numerical algorithm to approximate statistical solutions of multi-dimensional hyperbolic systems of conservation laws. To this end, we combined a high-resolution finite volume method (4.1) with a Monte Carlo sampling procedure to obtain Algorithm 1, that computes statistical solutions.

The task of proving convergence of the approximations $\mu_{t}^{{\Delta},M}$, generated by the Monte Carlo Algorithm, was rather intricate. First, we had to completely characterize an appropriate topology on the space of time-parameterized probability measures on $p$-integrable functions. This topology is based on the topology induced by the underlying correlation measures. We showed that this induced topology is equivalent to the weak topology on the space of probability measures on $L^{p}$, for any fixed time, but it also induces appropriate extensions when time is varied. The resulting compactness theorem 2.21 delineates the class of admissible observables that converge in this topology. Essentially, this theorem boils down to the convergence of time averages of (multi-point) statistical quantities of interest such as the mean, variance, multi-point correlation functions and structure functions (4.7). We believe that this topology on time-parameterized probability measures on $L^{p}$ and novel sufficient conditions for ensuring convergence in it, might have independent applications in probability theory and stochastic analysis.

Next, we proved in Theorem 4.1 that, under certain assumptions on the underlying finite volume schemes, the approximate statistical solutions converge, upon mesh refinement, in the aforementioned topology to a time-parameterized probability measure on $L^{p}$. A Lax–Wendroff theorem was proved, showing that the limit measure is indeed a statistical solution. Finally, a standard Monte Carlo convergence argument was used to guarantee convergence, under sample augmentation, of the approximations generated by the Monte Carlo Algorithm.

The assumptions in Theorem 4.1 include $L^{p}$ stability and the weak BV bound (4.5), which are satisfied by many existing high-resolution entropy stable finite volume schemes, such as the TeCNO schemes of [18]. On the other hand, we also required a subtle (approximate) scaling assumption (4.6) that is an analogue of the well-known scaling assumptions in Kolmogorov’s theories for homogeneous, isotropic turbulence [24]. Although this assumption was verified in all the numerical experiments, we were unable to prove it here. Hence, we have provided a conditional existence result for statistical solutions of multi-dimensional systems of conservation laws in this paper.

Moreover, we proposed an entropy condition and a notion of dissipative statistical solutions (Definition 3.2) and proved a weak-strong uniqueness result for these dissipative statistical solutions. This provides us with a conditional uniqueness result for statistical solutions, i.e., if they exist, then strong statistical solutions (Definition 3.4) are unique. In particular, we obtain short time existence and uniqueness results for dissipative statistical solutions.

We present extensive numerical experiments for the two-dimensional compressible Euler equations to illustrate our Monte Carlo Algorithm. The results validate the convergence analysis and demonstrate convergence, on mesh refinement and sample augmentation, for statistical quantities of interest such as the mean, variance, structure functions, one-point probability density functions and two-point joint probability density functions. Moreover, we observe convergence in appropriate Wasserstein distances (5.6) for the multi-point correlation measures. Summarizing the results of the numerical experiments, we conclude that one observes convergence of all interesting statistical observables in our framework. This should be contrasted to the state of the art, where deterministic quantities do not converge on mesh refinement [15]. Thus, we provide rigorous justification of the computability of statistical quantities of interest in the context of multi-dimensional systems of conservation laws.

Furthermore, we discover from the numerical experiments that

• •

The computed statistical solutions are remarkably stable with respect to different variations. In particular, we varied the amplitude of initial perturbations, the type of initial perturbations leading to different probability measures on $L^{p}$, and also the underlying numerical method. In all these cases, we observed that the computed statistical solutions were stable with respect to these perturbations. This observed stability augurs very well for identifying further constraints or admissibility criteria in order to obtain uniqueness of statistical solutions.

• •

The correlation structure of the statistical solutions seem to reach an equilibrium at significantly shorter time scales than the actual flow. This behavior is clearly seen in the variation of structure functions over time (Figure 18). Thus, statistical stationary might be reached much faster than the actual evolution would suggest.

• •

There seems to be a subtle gain in regularity, as measured by the decay exponent of the structure function, for the statistical solutions (Figure 18). Qualitatively, it seems that mixing of the underlying structures leads to a gain in regularity. This observed regularity needs to be studied further.

The results in this paper can be extended in different directions. On the theoretical side, a key question is whether the scaling assumption (4.6) can be proved for some numerical approximations or relaxed in an appropriate manner. This will pave the way for a unconditional global existence result for statistical solutions.

As formulated, (dissipative) statistical solutions are not necessarily unique. If we start with a deterministic initial data, i.e. setting $\bar{\mu}=\delta_{\bar{u}}$, then we can apply the construction of [11, 9] to obtain infinitely many entropy solutions corresponding to the initial data. Each deterministic solution defines a statistical solution. Thus, infinitely many statistical solutions are possible for the same (deterministic) initial data. On the other hand, numerical experiments strongly suggest that the computed statistical solutions are stable. Thus, we need to find further admissibility criteria to single out an unique statistical solution. Moreover, the observed gain in regularity might provide additional constraints to obtain uniqueness.

On the computational side, the Monte Carlo algorithm 1 can be very expensive, even prohibitively expensive in three space dimensions. Hence, it is imperative to consider alternatives to accelerate it. Alternatives such as the multi-level Monte Carlo method and quasi-Monte Carlo method are considered in a forthcoming publication. Another alternative would be to use deep learning, such as in [37], to accelerate the Monte Carlo algorithm.

Acknowledgements.

The research of SM was partially supported by the ERC consolidator grant ERC COG NN. 770880 COMANFLO. Numerical results in this paper were obtained from computations performed on the Swiss National supercomputing center (CSCS) Lugano, through production projects s665 and s839. FW was partially supported by NFR project no. 250302.

Appendix A Proof of Theorem 2.5

In order to characterize weak convergence of probability measures on $L^{p}(D;U)$, we will use a construction from [2, Section 5.1, pp. 106–107] which we summarize here. Assume that $X$ is Polish (i.e., a complete and separable metric space) and let $X_{0}\subset X$ be a dense, countable subspace. For fixed $q_{1},q_{2},q_{3}\in\mathbb{Q}$ with $q_{2},q_{3}\in(0,1)$ and every $v\in X_{0}$, define $H:X\to\mathbb{R}$ by

[TABLE]

The collection $\widehat{\Lambda}_{0}$ of all such “hat functions” $H$ is clearly countable. Let

[TABLE]

Then $\Lambda_{0}$ is also countable, and it can be shown that weak convergence of a sequence $(\mu_{n})_{n}$ is equivalent to

[TABLE]

(see the aforementioned reference).

We prove first Lemma 2.4.

Proof of Lemma 2.4.

If $g\in\mathcal{H}^{k,p}(D;U)$ and $u\in L^{p}(D;U)$ then

[TABLE]

If $g\in\mathcal{H}^{k,p}_{1}$ and $u,v\in L^{p}(D;U)$ then ($\hat{u}^{j}(x):=(u(x_{1}),\dots,u(x_{j-1}),u(x_{j+1}),\dots,u(x_{k}))$)

[TABLE]

whence $L_{g}$ is continuous (and in fact locally Lipschitz). ∎

Lemma A.1**.**

Let F\in\big{(}\mathcal{C}^{p}\big{)}^{k} for some $k\in\mathbb{N}$, i.e. $F(u)=(F_{1}(u),\dots,F_{k}(u))$ for $F_{1},\dots,F_{k}$ lying either in $\mathcal{C}^{p}$ or in $\mathcal{C}^{p}_{1}$. Let $P:\mathbb{R}^{k}\to\mathbb{R}$ be a polynomial in $k$ variables. Then $P\circ F$ lies in $\mathcal{C}^{p}$ or in $\mathcal{C}^{p}_{1}$, respectively.

Proof.

$\mathcal{C}^{p}$ is clearly closed under scalar multiplication, as $\alpha L_{g}=L_{\alpha g}$ for any $\alpha\in\mathbb{R}$. It is therefore enough to show that $\mathcal{C}^{p}$ ($\mathcal{C}^{p}_{1}$, resp.) is closed under addition and multiplication. Let $g_{i}:D^{k_{i}}\times U^{k_{i}}\to\mathbb{R}$ for $i=1,2$ be Carathéodory functions satisfying (2.3), (2.4). Then

[TABLE]

where $k=k_{1}+k_{2}$ and

[TABLE]

The function $g$ is readily seen to satisfy (2.3), and using (2.3) it can be seen that $g$ also satisfies (2.4) whenever $g_{1},g_{2}$ do so. This shows that $\mathcal{C}^{p}$ ($\mathcal{C}^{p}_{1}$, resp.) is closed under multiplication. To show that $L_{g_{1}}+L_{g_{2}}\in\mathcal{C}^{p}$ (or $\mathcal{C}^{p}_{1}$, respectively), assume that, say, $k_{1}\leqslant k_{2}$. If $k_{1}=k_{2}$ then clearly $L_{g_{1}}+L_{g_{2}}\in\mathcal{C}^{p}$ (or $\mathcal{C}^{p}_{1}$, respectively), so assume that $k_{1}<k_{2}$. Let $k_{0}=k_{2}-k_{1}$ and let $\varphi\in L^{1}(D^{k_{0}})$ satisfy $\varphi\geqslant 0$ and $\int_{D^{k_{0}}}\varphi(x)\,dx=1$. Define $g:D^{k_{2}}\times U^{k_{2}}\to\mathbb{R}$ by

[TABLE]

It is now straightforward to verify that $L_{g_{1}}+L_{g_{2}}=L_{g}$ and that $g$ satisfies (2.3), (2.4). ∎

Proof of Theorem 2.5.

We may assume that $B$ is closed. By Lemma 2.4, every $L_{g}\in\mathcal{C}^{p}_{1}$ is continuous and bounded on $B$. For any $L_{g}\in\mathcal{C}^{p}_{1}$ there is, by the Tietze extension theorem [22, Section 4.2], an $F\in C_{b}(L^{p}(D;U))$ satisfying $F=L_{g}$ on $B$, so we see that if $\mu_{n}\rightharpoonup\mu$ then \bigl{\langle}\mu_{n},L_{g}\bigr{\rangle}=\bigl{\langle}\mu_{n},F\bigr{\rangle}\to\bigl{\langle}\mu,F\bigr{\rangle}=\bigl{\langle}\mu,L_{g}\bigr{\rangle}.

Assume conversely that \bigl{\langle}\mu_{n},F\bigr{\rangle}\to\bigl{\langle}\mu,F\bigr{\rangle} for all $F\in\mathcal{C}^{p}_{1}$. We claim that the “hat functions” $H:L^{p}(D;U)\to\mathbb{R}$ defined in (A.1) can be approximated uniformly on $B$ by functions in $\mathcal{C}^{p}_{1}$. Let $v\in L^{p}(D;U)$ and $q_{1},q_{2},q_{3}\in\mathbb{Q}$ be as in (A.1), and let $R:=\sup_{u\in B}\|u-v\|_{L^{p}}^{p}<\infty$. Let \psi(s)=\min\big{(}q_{1}+q_{2}|s|^{1/p},\ q_{3}\big{)}, which is a continuous, bounded function on $\mathbb{R}$. Let $\varepsilon>0$ and let $P:\mathbb{R}\to\mathbb{R}$ be a polynomial such that $\|\psi-P\|_{C([-R,R])}<\varepsilon$. The function $u\mapsto\|u-v\|_{L^{p}}^{p}$ clearly lies in $\mathcal{C}^{p}_{1}$, so Lemma A.1 implies that also u\mapsto P\big{(}\|u-v\|_{L^{p}}^{p}\big{)} can be written as a function $L_{g}\in\mathcal{C}^{p}_{1}$. For any $u\in B$ we then have

[TABLE]

This proves the claim.

Next, let $F\in\Lambda_{0}$, where $\Lambda_{0}$ is given by (A.2). Let $\varepsilon>0$, let $g_{1},\dots,g_{m}$ be such that $\|H_{i}-L_{g_{i}}\|_{C_{b}(B)}<\varepsilon$ for $i=1,\dots,m$, and let $P:\mathbb{R}^{m}\to\mathbb{R}$ be a polynomial such that

[TABLE]

By Lemma A.1 we can write P\circ\big{(}L_{g_{1}},\dots,L_{g_{m}}\big{)}=L_{g} for some $L_{g}\in\mathcal{C}^{p}_{1}$. We conclude that

[TABLE]

by the 1-Lipschitz continuity of the min-function and the approximation properties of $g_{1},\dots,g_{m}$ and $P$. We can conclude that (A.3) holds, and hence $\mu_{n}\rightharpoonup\mu$ weakly. ∎

Appendix B The compactness theorem

Throughout this appendix we will use the cutoff functions

[TABLE]

(for some $R\geqslant 1$), defined for $s\in\mathbb{R}$ and $v\in U$, where $U\subset\mathbb{R}^{N}$ as before. Note that $\zeta_{R}(v)=v$ for $|v|<R$ and $\zeta_{R}(v)=0$ for $|v|>2R$, and that $\|\zeta_{R}\|_{\operatorname{Lip}}\leqslant 2$.

For a function $u\in L^{p}(D)$ we define the modulus of continuity

[TABLE]

Lemma B.1**.**

If $u\in L^{q}(D,U)$ for some $q\in[1,\infty)$ then

[TABLE]

whenever $R\geqslant\frac{1}{|B_{r}|^{1/q}}\|u\|_{L^{q}(D)}$.

Proof.

Denote u_{r}(x):=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.31105pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.50664pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.4333pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-3.96664pt}}\!\int_{B_{r}(0)}u(x+z)\,dz. Then $|u_{r}(x)|\leqslant\frac{1}{|B_{r}|^{1/q}}\|u\|_{L^{q}(D)}$ for every $x\in D$, so that if $R\geqslant\frac{1}{|B_{r}|^{1/q}}\|u\|_{L^{q}(D)}$ then $u_{r}=\zeta_{R}\circ u_{r}$. It follows that

[TABLE]

where we have used the fact that $\|u_{r}-u\|_{L^{q}}\leqslant\omega_{r}^{q}(u)^{1/q}$. ∎

Proof of Theorem 2.13.

For every $k\in\mathbb{N}$, the sequence $(\nu^{k}_{n})_{n\in\mathbb{N}}\subset\mathcal{H}^{k*}_{0}(D;U)$ is bounded with norm $\|\nu^{k}_{n}\|_{\mathcal{H}^{k*}}\equiv 1<\infty$, and so has a weak*-convergent subsequence. Thus, we can extract a diagonal subsequence $(n_{j})_{j\in\mathbb{N}}$ such that $\nu^{k}_{n_{j}}\overset{*}{\rightharpoonup}\nu^{k}\in\mathcal{H}^{k*}_{0}(D;U)$ for every $k\in\mathbb{N}$. For the sake of notational simplicity we denote $n=n_{j}$ for the remainder of this proof.

We start by showing (iv). Let $\kappa$ and $\varphi$ be nonnegative functions, as prescribed. We may assume that \liminf_{n\to\infty}\bigl{\langle}\nu^{k}_{n},g\bigr{\rangle}_{\mathcal{H}^{k}}<\infty, for otherwise there is nothing to prove. For $R>0$, let $\kappa_{R}\in C_{c}(U^{k})$ and $\varphi_{R}\in L^{1}(D^{k})$ be the functions $\kappa_{R}(\xi)=\kappa(\xi)\theta_{R}(\xi)$ and $\varphi_{R}(x)=\varphi(x)\theta_{R}(x)$. Then $\kappa_{R}\to\kappa$ and $\varphi_{R}\to\varphi$ pointwise almost everywhere as $R\to\infty$. If $g_{R}(x,\xi):=\varphi_{R}(x)\kappa_{R}(\xi)$ then clearly $g_{R}\to g$ as $R\to\infty$ almost everywhere, and $g_{R}\in\mathcal{H}^{k}_{0}(D;U)$ for every $R$. We then find that

[TABLE]

which is (2.13).

Since $\nu^{k}_{n;x}\geqslant 0$ for all $k,n$ and a.e. $x\in D^{k}$, we have $\nu^{k}_{x}\geqslant 0$ for all $k$ and a.e. $x\in D^{k}$. To see that $\nu^{k}_{x}(U^{k})=1$ for a.e. $x\in D^{k}$, let $A_{1},\dots,A_{k}\subset D$ be bounded Borel sets of positive Lebesgue measure and define $A=A_{1}\times\dots\times A_{k}$. Then since $\nu^{k}_{n;x}$ are probability measures and converge weak* as we have just shown,

[TABLE]

Conversely,

[TABLE]

where we have first used Chebyshev’s inequality and then the uniform $L^{p}$ bound (2.11). Passing to $R\to\infty$ in the two estimates above shows that

[TABLE]

Since $A_{1},\dots,A_{k}$ were arbitrary, it follows that $\nu^{k}_{x}(U^{k})=1$ for a.e. $x\in D^{k}$.

The limit $\bm{\nu}$ clearly satisfies the symmetry condition. For consistency, we need to show that for any $f\in C_{0}(U^{k-1})$ and $\varphi\in C_{c}(D^{k})$, we have

[TABLE]

Define $\bar{f}_{R}\in C_{0}(U^{k})$ by $\bar{f}_{R}(\xi_{1},\dots,\xi_{k})=f(\xi_{1},\dots,\xi_{k-1})\theta_{R}(\xi_{k}).$ We estimate the difference between the left- and right-hand sides of (B.2) by $E_{1}+E_{2}+E_{3}+E_{4}$, where

[TABLE]

where $\hat{x}_{k}=(x_{1},\dots,x_{k-1})$. The fact that $E_{1}\to 0$ as $R\to\infty$ follows from the Dominated Convergence Theorem. For every $R>0$ we have $E_{2}$, $E_{4}\to 0$ as $n\to\infty$ by the weak* convergence of $\nu^{k}_{n}$ and $\nu^{k-1}_{n}$, respectively. For $E_{3}$ we can write

[TABLE]

as $R\to\infty$ (as in the proof of (B.1)). Thus, for a given $\varepsilon>0$ we can first select $R$ such that $E_{1},E_{3}<\varepsilon$, and then select $n$ such that $E_{2},E_{4}<\varepsilon$. This proves (B.2).

The $L^{p}$ bound (ii) follows from (2.13) with $k=1$, $\varphi\equiv 1$ and $\kappa(\xi)=|\xi|^{p}$, since \bigl{\langle}\nu^{1},|\xi|^{p}\bigr{\rangle}_{\mathcal{H}^{1}}\leqslant\liminf_{n}\bigl{\langle}\nu^{1}_{n},|\xi|^{p}\bigr{\rangle}_{\mathcal{H}^{1}}\leqslant c^{p}, by (2.11). To prove the diagonal continuity property (iii) of $\bm{\nu}$ we use (2.13) with $k=2$, $\varphi(x,y)=\frac{\mathbf{1}_{B_{r}(x)}(y)}{|B_{r}(x)|}$ and $\kappa(\xi_{1},\xi_{2})=|\xi_{1}-\xi_{2}|^{p}$, and note that \omega_{r}^{p}(\nu^{2})=\bigl{\langle}\nu^{2},g\bigr{\rangle}_{\mathcal{H}^{2}}. The diagonal continuity of $\nu^{2}$ then follows from $\omega_{r}^{p}(\nu^{2})\leqslant\liminf_{n\to\infty}\omega_{r}^{p}(\nu^{2}_{n})\to 0$ as $r\to 0$, by (2.12).

We prove (v) for $k=1$ and give a sketch of the general case. Let $g\in\mathcal{H}^{1,p}_{1}(D;U)$, that is, $g$ is a Carathéodory function such that $|g(x,\xi)|\leqslant\varphi_{1}(x)+\varphi_{0}|\xi|^{p}$ for some $\varphi_{0}>0$ and nonnegative $\varphi_{1}\in L^{1}(D)$, as well as $|g(x,\xi)-g(y,\zeta)|\leqslant C\max(|\xi|,|\zeta|)^{p-1}|\xi-\zeta|$ for some $C>0$. Define

[TABLE]

The sequence $\{u_{n}\}_{n\in\mathbb{N}}$ satisfies

[TABLE]

so together with the compactness of $D$ and the uniform bound $\|u_{n}\|_{L^{1}(D)}\leqslant\|\varphi_{1}\|_{L^{1}(D)}+\varphi_{0}c^{p}$ we can apply Kolmogorov’s Compactness Theorem [29, Theorem A.8], [43] to conclude that $\{u_{n}\}_{n\in\mathbb{N}}$ is precompact in $L^{1}(D)$. Hence, there exists a subsequence $\{u_{n_{l}}\}_{l\in\mathbb{N}}$ and some $\bar{u}\in L^{1}(D)$ such that $u_{n_{l}}\to\bar{u}$ as $l\to\infty$ in $L^{1}(D)$. From the weak* convergence of $\nu^{1}_{n}$ we know that $u_{n}^{R}\rightharpoonup u^{R}$ as $n\to\infty$ weakly in $L^{1}(D)$ for every $R>0$. Lebesgue’s dominated convergence theorem implies that $u^{R}\to u$ as $R\to\infty$ in $L^{1}(D)$.

We claim that $u^{R}_{n}\to u_{n}$ as $R\to\infty$ in $L^{1}(D)$, uniformly in $n$. Indeed, choosing $R>0$ such that $R\geqslant\frac{1}{|B_{r}|^{1/p}}M$ (where $M$ is the constant in (2.14)), we get

[TABLE]

The first term can be bounded by

[TABLE]

where we in the second inequality used Chebyshev’s inequality and in the third inequality the uniform $L^{p}$ bound (2.11). The above converges to 0 as $R\to\infty$, uniformly in $n$. For the second term we have

[TABLE]

where the first inequality follows from Lemma B.1 with $q=1$, the second inequality follows an estimate similar to (B.3), and the third inequality is Hölder’s inequality. The final term above vanishes as $R\to\infty$ uniformly in $n$, by the uniform diagonal continuity assumption (2.12). It follows that for any $\psi\in L^{\infty}(D)$,

[TABLE]

all of which vanish as $R\to\infty$ and $l\to\infty$. We conclude that $\bar{u}=u$, whence $u_{n_{l}}\to u$ as $l\to\infty$. By the uniqueness of the limit $u$, we get convergence of the whole sequence: $u_{n}\to u$ in $L^{1}(D)$ as $n\to\infty$.

For general $k\in\mathbb{N}$ we prove only that u_{n}(x):=\bigl{\langle}\nu^{k}_{n;x},g(x,\xi)\bigr{\rangle} satisfies a bound of the form (B.3), and leave the rest to the reader. We write first

[TABLE]

(cf. Definition 2.2 with $\hat{\xi}^{l}=(\xi_{1},\dots,\xi_{l-1},\xi_{l+1},\dots,\xi_{k})$). Consider, say, the last summand above:

[TABLE]

where the second-last inequality follows from Hölder’s inequality and the boundedness of $L_{h}$ (cf. the proof of Lemma 2.4), and the last inequality follows from (2.14). By the uniform diagonal continuity of $\bm{\nu}_{n}$, the above vanishes uniformly as $r\to 0$. The rest of the proof follows as in the proof for $\nu^{1}$. ∎

Appendix C Time-dependent correlation measures

Proof of Lemma 2.17.

By assumption, (t,x)\mapsto\bigl{\langle}\nu^{k}_{t,x},g\bigr{\rangle} is measurable for all $g\in C_{0}(U^{k})$, so by separability of $\mathcal{H}^{k}_{0}(D;U)$, the map t\mapsto\int_{D^{k}}\bigl{\langle}\nu^{k}_{t,x},g(x)\bigr{\rangle}\,dx is also measurable for any $g\in\mathcal{H}^{k}_{0}(D;U)$. Let $E_{0}\subset\mathcal{H}^{k}_{0}(D;U)$ be a countable, dense subset of the unit sphere in $\mathcal{H}^{k}_{0}(D;U)$, and let

[TABLE]

a set whose complement $[0,T)\setminus\mathcal{T}$ has zero Lebesgue measure. The set $E:=\textrm{span}\,E_{0}$ is dense in $\mathcal{H}^{k}_{0}(D;U)$, and every $t\in\mathcal{T}$ is still a Lebesgue point for t\mapsto\int_{D^{k}}\bigl{\langle}\nu^{k}_{t,x},g(x)\bigr{\rangle}\,dx whenever $g\in E$. For every $t\in\mathcal{T}$ we define the functional $\rho(t):E\to\mathbb{R}$ by

[TABLE]

Then $\rho(t)$ is linear: For any $g,h\in E$ and $\alpha\in\mathbb{R}$ we have

[TABLE]

Moreover, $\rho(t)$ is continuous, as

[TABLE]

It follows that for every $t\in\mathcal{T}$ the functional $\rho(t)$ can be extended uniquely to a continuous linear functional on $\mathcal{H}^{k}_{0}(D;U)$. The remaining claims in the lemma are now readily checked. ∎

Appendix D Proof of Theorem 4.6

We write

[TABLE]

where $\mathcal{S}^{{\Delta}}$ is the numerical evolution operator; see Section 4.1.

Proof.

Let $\mu^{{\Delta}}_{t}$ be the approximate statistical solution, defined in (4.3). Let $k\in\mathbb{N}$, and let $\varphi\in C_{c}^{\infty}(D^{k})$ be a test function. Fix ${\mathbf{i}}=({\mathbf{i}}_{1},\ldots,{\mathbf{i}}_{k})\in(\mathbb{Z}^{d})^{k}$, By changing the order of integration, we have

[TABLE]

where

[TABLE]

Since $\{\bar{u}^{\cdot,{\Delta}}_{{\mathbf{j}}}\}_{{\mathbf{j}}\in\mathbb{Z}^{d}}$ solves (4.1) with initial data $\left\{\bar{u}_{{\mathbf{j}}}\right\}_{{\mathbf{j}}\in\mathbb{Z}^{d}}$ for every $\bar{u}\in L^{1}(D,U)$, we get that

[TABLE]

Multiplying by ${\Delta}^{kd}$ and summing over $(\mathbb{Z}^{d})^{k}$, we get

[TABLE]

For $1\leqslant n\leqslant k$ and $1\leqslant m\leqslant d$, summation by parts gives

[TABLE]

We furthermore have that

[TABLE]

and by the Lipschitz continuity (4.2), we have

[TABLE]

By Hölder’s inequality we get

[TABLE]

By the weak BV assumption (4.18), we get

[TABLE]

We furthermore note that since $\varphi\in C_{c}^{\infty}$, we have

[TABLE]

Hence we get

[TABLE]

which completes the proof. ∎

Appendix E Proof of (5.7)

Proposition E.1**.**

Let $D\subset\mathbb{R}^{d}$ be bounded and let $\mu,\tilde{\mu}\in\mathcal{P}(L^{p}(D;U))$ both have bounded $p$th moment. Then

[TABLE]

for some $C=C(p,k,|D|)$, where $(\nu^{k})_{k\in\mathbb{N}},(\tilde{\nu}^{k})_{k\in\mathbb{N}}$ are the correlation measures corresponding to $\mu,\tilde{\mu}$.

Proof.

Let $\pi\in\mathcal{P}(L^{p}(D;U)^{2})$ be an optimal transport plan for $W_{p}(\mu,\tilde{\mu})$, and let $(\sigma^{k})_{k\in\mathbb{N}}\in\mathcal{L}^{p}(D;U^{2})$ be its corresponding correlation measure. We claim that the marginals of $\sigma^{k}_{x}$ are $\nu^{k}_{x}$ and $\tilde{\nu}^{k}_{x}$, respectively, for a.e. $x\in D^{k}$. Indeed, if $g\in\mathcal{H}^{k}_{0}(D;U)$ then, denoting $\bar{g}(x,\xi,\zeta):=g(x,\xi)$,

[TABLE]

A similar computation holds for the second marginal. Since $g$ is arbitrary, the claim follows. We can therefore estimate

[TABLE]

∎

Appendix F Computing the Wasserstein distance for sums of Diracs

In Section 5, we computed the Wasserstein distance between different numerical solutions. The numerical solutions are always sums of Diracs, which simplifies the computations greatly.

For two $k$-th correlation marginals $\nu^{k,1}$ and $\nu^{k,2}$, we are interested in the distance

[TABLE]

F.1. Computing the Wasserstein distance between the first correlation marginals

To compute the first correlation marginals, we use the function

scipy.stats.wasserstein_distance in the scipy module for Python [30] to compute the Wasserstein distance. The function computes the Wasserstein distance by going through the CDF, consult the scipy documentation for more information.

We approximate the spatial integral as a sum over all the volume averages.

F.2. Computing the Wasserstein distance between the second correlation marginals

For the second correlation marginals, we use the function ot.emd in the POT module for Python [20]. This function uses the Hungarian algorithm [33] to compute the Wasserstein distance between sums of Diracs.

The spatial integral is approximated by taking 10 spatial points in each direction, in other words, we use $10\;000$ evaluation of the Wasserstein distance.