\section{Introduction}

Conservation laws are widely used to model problems in various applications, such as fluid and gas dynamics \cite{Andreianov2014Approximating, Jaffre1996Numerical}, hydrodynamics \cite{Chen2008Hyperbolic}, nonlinear elasticity, semiconductor simulation and many more. We refer to the books of Dafermos \cite{Dafermos2000Hyperbolic} and LeVeque \cite{LeVeque1992Numerical} and the references therein for more applications.
There are two natural extensions of models based on conservation laws:
\begin{enumerate}[(i)]
	\item One may want to allow the flux function to depend (possibly discontinuously) on the state space, e.g., to model heterogeneities or fractures in a porous medium for subsurface or two-phase flow simulations \cite{Andreianov2014Approximating, Jaffre1996Numerical, Gimse1992Solution}.
	This approach is also used for many other models in a variety of fields, e.g.,  vehicular traffic flows \cite{Andreianov2010Finite, Coclite2005Traffic}, statistical mechanics \cite{Chen2008Hyperbolic}, sedimentation \cite{Buerger2004Well} and kidney physiology \cite{Perthame2014simple}.
	\item One might be interested in including stochastic terms, e.g., a stochastic flux function \cite{Holden1997Conservation, Lions2013Scalar, Lions2014Scalar} or a stochastic source \cite{Karlsen2017stochastic, Kim2003stochastic}, to account for uncertainties or insufficient data in the underlying model.
\end{enumerate}
To combine both extensions, we consider a stochastic conservation law, where the flux function depends discontinuously on the state space.
In this paper, we are concerned with well-posedness of the solution to such a stochastic conservation law.
A special case of the considered stochastic conservation law are flux functions, in which a discontinuous random coefficient is incorporated.
This random coefficient is constructed via a continuous part and a discontinuous part, which is inspired by the unique characterization of Lévy processes via the Lévy-Khintchine formula \cite{Applebaum2009Levy, Sato1999Levy}. By this formula, every Lévy process can be seen as the composition of three independent components, namely, a drift term, a Brownian motion and a pure jump process.

The paper is structured as follows:
In the remainder of this section, we discuss existing results on conservation laws with discontinuous flux function in Section \ref{ssec:intro_discontinuous_CLs}, before stating existing results on stochastic conservation laws in Section \ref{ssec:intro_SCLs}.
In the following section, we state the considered problem and show our main result: the well-posedness of pathwise solutions.
In Section \ref{sec:stochJumpCoeff}, we define a random coefficient and its numerical approximation.
Afterwards, we introduce a numerical discretization of the stochastic conservation law with discontinuous random coefficient and investigate the convergence of the approximations in numerical experiments.

\subsection{Conservation laws with spatial discontinuities}
\label{ssec:intro_discontinuous_CLs}

We are interested in scalar conservation laws with a flux that depends discontinuously on the spatial variable.
For $t \in \mathbb{T}$, where $\mathbb{T} \subset \mathbb{R}$ denotes the time interval, the conservation law is given by
\begin{equation}
	\label{eq:deterministicConservationLaw}
	\begin{aligned}
		u_{t} + \mathfrak{\solFlux}(x, u)_{x} &= 0 &&\text{in } \mathbb{R} \times \mathbb{T} \ , \\
		u(x, 0) &= \solution_{0} (x) &&\text{in } \mathbb{R} \ ,
	\end{aligned}
\end{equation}
where $\solution_{0} (x) \in L^{\infty}(\mathbb{R})$ denotes the initial condition of the partial differential equation.
Furthermore, the flux $\mathfrak{\solFlux}$ is assumed to depend discontinuously on space.
In case the flux function $\mathfrak{\solFlux}$ is sufficiently smooth, well-posedness of \eqref{eq:deterministicConservationLaw} is well established \cite{Bouchut1998Krukovs, Kruzkov1970First}.
However, even with smooth initial data and a sufficiently smooth flux function, the solution might develop discontinuities.
Therefore, weak solutions are sought, but they are not unique without an additional admissibility criterion selecting the physical solution.
In case of a smooth space dependency of the flux $\mathfrak{\solFlux}$ the entropy condition was described by Kru\v{z}kov \cite{Kruzkov1970First}.
A function $u \in L^{\infty}(\mathbb{R} \times \mathbb{T}; \mathbb{R})$ is an entropy solution of \eqref{eq:deterministicConservationLaw} if for all $m \in \mathbb{R}$ it satisfies the following inequality in the distributional sense:
\begin{equation}
	\label{eq:KruzkovEntropyCondition}
	\frac{\partial}{\partial \timeVar} \abs{u - m} + \frac{\partial}{\partial \spaceVar} \left( \sign(u - m) (\mathfrak{\solFlux}(x, u) - \mathfrak{\solFlux}(x, m)) \right) + \sign(u - m) \frac{\partial}{\partial \spaceVar} \mathfrak{\solFlux}(x, m) \leq 0 \ .
\end{equation}
For discontinuous flux functions, the Kru\v{z}kov entropy condition \eqref{eq:KruzkovEntropyCondition} does not make sense due to the term $\sign(u - m) \frac{\partial}{\partial \spaceVar} \mathfrak{\solFlux}(x, m)$.
To overcome this difficulty, a variety of criteria were introduced to select the appropriate weak solution.
However, for some problems, different selection criteria might lead to different solutions.
Therefore, the solution we are seeking might depend on the context.

A common approach to define a selection criterion is to require the usual entropy condition, when the flux is continuous, and impose additional constraints at the discontinuities (see, e.g., \cite{Adimurthi2005Godunov, Bulicek2011scalar, Ghoshal2011Existence}).
A natural choice of such a condition is the continuity of the flux at the heterogeneities,  given by the Rankine-Hugoniot condition
\begin{equation}
	\label{eq:Rankine-Hugoniot}
	\mathfrak{\solFlux}(x^-, u(x^-, t)) = \mathfrak{\solFlux}(x^+, u(x^+, t)) \quad t \in \mathbb{T} .
\end{equation}
However, this condition is not sufficient to select an unique solution.
Note that for discontinuous flux functions this procedure requires the existence of traces along the discontinuities, which implies the spatial points to occur at isolated points.
For an exhaustive review of different selection criteria, the reader is referred to \cite{Andreianov2011theory}, where each existing admissibility criterion is associated to a so-called \textit{Germ}.

\subsubsection{Audusse-Perthame adapted entropy formulation}

The theory of adapted-entropy solutions for conservation laws was introduced by Audusse and Perthame \cite{Audusse2005Uniqueness} in 2005, based on the ideas of Baiti and Jenssen \cite{Baiti1997Well} from 1997.
The provided framework is formulated, such that no additional interface conditions are imposed.
As a result, the existence of traces is not required, since the discontinuity points are not distinguished separately. This makes it possible to include an infinite number of discontinuity points in the flux function.

The main idea of the framework is to adapt the classical Kru\v{z}khov entropy condition \eqref{eq:KruzkovEntropyCondition} in such a way that the troublesome term vanishes.
Therefore, instead of the classical Kru\v{z}khov \cite{Kruzkov1970First} entropies $(\abs{u - m})_{m \in \mathbb{R}}$, adapted entropies $(\abs{u - \steadyStateSolution{}})_{\alpha \in \mathbb{R}}$ are considered.
Here, $\steadyStateSolution{}$ denotes the unique solution of the problem
\begin{equation}
	\label{eq:steadyStateFlux}
	\mathfrak{\solFlux}(x, \steadyStateSolution{}) = \alpha \quad \text{for a.e. }x \in \mathbb{R} \ .
\end{equation}
The existence and uniqueness of $\steadyStateSolution{}$ in \eqref{eq:steadyStateFlux} is guaranteed by the following assumptions:
\begin{assumption}[Assumption on the flux function: adapted entropy formulation]
	\label{as:deterministicFlux}
	\begin{enumerate}[label=(A-\arabic*)]
		\item \label{as:A-1} The flux function $\mathfrak{\solFlux}(x, \rho)$ is continuous at all points of $(\mathbb{R} \setminus \mathcal{N}) \times \mathbb{R}$, where $\mathcal{N} \subset \mathbb{R}$ is a closed set of measure zero that contains the discontinuity points of $\mathfrak{\solFlux}(\cdot, \rho)$.
		\item \label{as:A-2} There exist two functions $g_{-}, g_{+} \in C^{0}(\mathbb{R})$ such that for all $x \in \mathbb{R}$ it holds that $g_{-} (\rho) \leq |\mathfrak{\solFlux}(x, \rho)| \leq g_{+} (\rho)$, where $g_{-}$ is a non-negative (non-strictly) decreasing then increasing function which satisfies $|g_{-}(\pm \infty)| = +\infty$.
		\item \label{as:A-3} There exists a function $\solution_M (x): \mathbb{R} \rightarrow \mathbb{R}$ and a constant $M_0$ such that for $x \in \mathbb{R} \setminus \mathcal{N}$, $\mathfrak{\solFlux}(x, \cdot)$ is a locally Lip\-schitz one-to-one function from $(-\infty, \solution_M(x)]$ and $[\solution_M(x), +\infty)$ to $[M_0, +\infty)$ (or $(-\infty, M_0]$) that satisfies $\mathfrak{\solFlux}( x, \solution_M(x)) = 0$ and with common Lipschitz constant $L_I$ for all $x \in \mathbb{R} \setminus \mathcal{N}$ and all $u \in I$, where $I \subset \mathbb{R}$ is any bounded interval.
	\end{enumerate}
	Alternatively, instead of \ref{as:A-3} we may consider the assumption
	\begin{enumerate}
		\item[\mylabel{as:A-3alt}{(A-3')}] For $x \in \mathbb{R} \setminus \mathcal{N}$, the flux function $\mathfrak{\solFlux} (x, \cdot)$ is a locally Lipschitz one-to-one function from $\mathbb{R}$ to $\mathbb{R}$ with common Lipschitz constant $L_I$ for all $x \in \mathbb{R} \setminus \mathcal{N}$ and all $u \in I$, where $I \subset \mathbb{R}$ is any bounded interval.
	\end{enumerate}
\end{assumption}

In case the flux function satisfies Assumptions \ref{as:A-1}\--\ref{as:A-3}, for any constant $\alpha \in [M_0, \infty)$ (or $\alpha \in (-\infty, M_0]$), there exist two steady-state solutions $\steadyStateSolution{+}$ from $\mathbb{R}$ to $[\solution_M (x), \infty)$ and $\steadyStateSolution{-}$ from $\mathbb{R}$ to $(-\infty, \solution_M (x)]$ of \eqref{eq:deterministicConservationLaw} satisfying \eqref{eq:steadyStateFlux}.
In case of Assumptions \ref{as:A-1}\--\ref{as:A-3alt}, the two steady state solutions coincide, which is even simpler.

Now, the adapted entropies $(\abs{u - \steadyStateSolution{}})_{\alpha \in \mathbb{R}}$ allow us to define adapted entropy solutions.
\begin{definition}[adapted entropy solution]
	\label{def:adaptedEntropySolution}
	A function $u \in L^{\infty}(\mathbb{R} \times \mathbb{T}) \cap C^{0}(\mathbb{T}, L^{1}_{\text{loc}}(\mathbb{R}))$ is an adapted entropy solution of \eqref{eq:deterministicConservationLaw} on $\mathbb{R} \times \mathbb{T}$, provided that for each $\alpha \in [M_0, \infty)$ (or $\alpha \in (-\infty, M_0]$) and the corresponding two steady-state solutions $\entropyConstant^{\pm}_{\steadyStatePar}(\spaceVar)$ of \eqref{eq:deterministicConservationLaw}, the following inequality holds in the sense of distributions:
	\begin{equation}
		\label{eq:adaptedEntropyCondition}
		\frac{\partial}{\partial \timeVar} \abs{u(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\spaceVar)} + \frac{\partial}{\partial \spaceVar} \left( \sign \left( u(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\spaceVar) \right) (\mathfrak{\solFlux}(x, u(x, t)) - \mathfrak{\solFlux}(x, \entropyConstant^{\pm}_{\steadyStatePar}(\spaceVar))) \right) \leq 0 \ .
	\end{equation}
\end{definition}

This formulation guarantees the uniqueness of solutions, since the \textit{doubling of variables} approach can be applied.
The existence was shown in \cite{Chen2008Hyperbolic} via the reduction of measure-valued solutions.
However, to the best of the author's knowledge, extending the adapted entropy formulation to multidimensional conservation laws is still an open problem.

\subsection{Stochastic conservation laws}
\label{ssec:intro_SCLs}

The field of stochastic conservation laws has attracted considerable attention in the past 20 years.
The concept of stochastic weak and entropy solutions was initially introduced by Holden and Risebro in \cite{Holden1997Conservation}, who considered nonlinear hyperbolic problems with time-dependent stochastic source terms.
This concept was further developed by Kim in \cite{Kim2003stochastic} and extended to conservation laws with multiplicative noise, see \cite{Feng2008Stochastic} for Gaussian and \cite{Biswas2015Conservation} for L\'{e}vy noise.
In \cite{Bauzet2012Cauchy}, Bauzet et al. used the concept of Kru\v{z}kov's entropy formulation and measure-valued solutions to establish existence and uniqueness of entropy solutions for multi-dimensional nonlinear conservation laws with multiplicative stochastic perturbation.
Recently, Karlsen and Storr{\o}sten generalized this approach by allowing the stochastic Kru\v{z}kov entropy condition to contain Malliavin differentiable random variables \cite{Karlsen2017stochastic}.
They also provided existence and uniqueness results for this extended framework.

In 2013 Lions \cite{Lions2013Scalar} considered scalar conservation laws with rough stochastic fluxes and extended this in \cite{Lions2014Scalar} to the case, where the rough stochastic flux function is space dependent.
The field of scalar conservation laws with rough fluxes combined with stochastic forcing was also considered by Hofmanov\'{a} in \cite{Hofmanova2016Scalar}, who proved well-posedness for the corresponding kinetic formulation and its solution.

Stochastic systems of conservation laws were considered by Mishra and Schwab in \cite{Mishra2012Sparse}, who introduced the uncertainty via random initial data.
In \cite{Poette2009Uncertainty} Po{\"e}tte et al. considered hyperbolic systems with a randomly perturbed flux function and showed the hyperbolicity of a stochastic Galerkin representation.
The stochastic Galerkin representation was later extended by \cite{Tryoen2012Adaptive} and \cite{Buerger2014hybrid}.
A comprehensive discussion on current numerical methods for hyperbolic systems of conservation laws can be found in \cite{Abgrall2017Uncertainty}.

In \cite{Mishra2016Numerical}, Mishra et al. considered the numerical approximation of scalar conservation laws in several space dimensions with random flux functions and also showed convergence analysis for the multilevel Monte Carlo finite volume method \cite{Mishra2012Sparse, Mishra2012Multi, Mishra2012Multilevel}.
Another numerical approximation was introduced by Risebro et. al. \cite{Risebro2016Multilevel}, who proposed a multilevel Monte Carlo front tracking algorithm for stochastic scalar conservation laws with bounded random flux function.
This approach was generalized in \cite{Koley2017multilevel} to the case of random degenerate convection-diffusion equations.

Recently, Barth and Stein \cite{Barth2019stochastic} considered semilinear hyperbolic stochastic partial differential equations, which are driven by L\'{e}vy noise and provided a fully discrete numerical approximation for this problem.

\section{Stochastic scalar conservation laws with discontinuous flux functions}

Let $(\stochDomain, \sigmaAlgebra, \mathbb{P})$ be a complete probability space and let $\mathbb{T} := [0, T] \subset \mathbb{R}_{\geq 0}$ be a time interval for some $T \in \mathbb{R}_{>0}$.
We consider a scalar, one-dimensional conservation law of the form
\begin{equation}
	\label{eq:stochasticConservationLaw}
	\begin{aligned}
		u_{t} + \mathfrak{\solFlux}(\omega, x, u)_{x} &= 0 &&\text{in } \Omega \times \mathbb{R} \times \mathbb{T} \ , \\
		u(\omega, x, 0) &= \solution_{0} (\omega, x) &&\text{in } \Omega \times \mathbb{R} \ ,
	\end{aligned}
\end{equation}
where $\solution_{0} \in L^p(\Omega, L^{\infty}(\mathbb{R}))$ denotes the stochastic initial condition of the partial differential equation.
We impose the following assumptions on the stochastic flux function:

\begin{assumption}[Pathwise assumption on the stochastic flux function]
	\label{as:stochFlux}
	We assume the stochastic discontinuous flux function to satisfy the following assumptions for every $\omega \in \Omega$:
	\begin{enumerate}[label=(B-\arabic*)]
		\item \label{as:B-1} The flux function $\mathfrak{\solFlux}(\omega, \cdot, \cdot)$ is continuous at all points of $(\mathbb{R} \setminus \mathcal{N}(\omega)) \times \mathbb{R}$ with $\mathcal{N}(\omega) \subset \mathbb{R}$ a closed set of measure zero that contains the discontinuity points of $\mathfrak{\solFlux} (\omega, \cdot, \rho)$.
		\item \label{as:B-2} There exist two functions $g_{-} (\omega, \cdot), g_{+} (\omega, \cdot) \in C^{0}(\mathbb{R})$ such that for all $x \in \mathbb{R}$ it holds that $g_{-} (\omega, \rho) \leq |\mathfrak{\solFlux}(\omega, x, \rho)| \leq g_{+} (\omega, \rho)$, where $g_{-} (\omega, \cdot)$ is a decreasing then increasing function with $g_{-}(\omega, \rho) \geq 0$ and $|g_{-}(\omega, \pm \infty)| = +\infty$.
		\item \label{as:B-3} There exists a function $\solution_M (\omega, x): \Omega \times \mathbb{R} \rightarrow \mathbb{R}$ and a constant $M_0 \in \mathbb{R}$ such that for  $x \in \mathbb{R} \setminus \mathcal{N}(\omega)$, the flux function $\mathfrak{\solFlux}(\omega, x, \cdot)$ is locally Lipschitz and one-to-one from $(-\infty, \solution_M (\omega, x)]$ and $[\solution_M(\omega, x), +\infty)$ to $[M_0, +\infty)$ (or $(-\infty, M_0]$). Further it satisfies $\mathfrak{\solFlux}(\omega, x, \solution_M(\omega, x)) = M_0$ with common Lipschitz constant $L_{\lipschitzInterval}$ for all $x \in \mathbb{R} \setminus \mathcal{N}(\omega)$ and all $u \in I$, where $I \subset \mathbb{R}$ is any bounded interval.
	\end{enumerate}

	Alternatively, instead of \ref{as:B-3} we may consider the assumption
	\begin{enumerate}
		\item[\mylabel{as:B-3alt}{(B-3')}] For $x \in \mathbb{R} \setminus \mathcal{N} (\omega)$, the flux function $\mathfrak{\solFlux} (\omega, x, \cdot)$ is locally Lipschitz and one-to-one from $\mathbb{R}$ to $\mathbb{R}$ with common Lipschitz constant $L_{\lipschitzInterval}$ for all $x \in \mathbb{R} \setminus \mathcal{N}(\omega)$ and all $u \in I$, where $I \subset \mathbb{R}$ is any bounded interval.
	\end{enumerate}
\end{assumption}

\begin{example}[Multiplicative flux]
	\label{ex:multiplicativeFlux}

	One example of stochastic flux functions satisfying Assumption \ref{as:stochFlux} is
	\begin{equation}
		\label{eq:multiplicativeFlux}
		\mathfrak{\solFlux}(\omega, x, \rho) = a (\omega, x) f (\rho) .
	\end{equation}
	Here, $a(\omega, x)$ is continuous in $x \in \mathbb{R}$ except on a closed set $\mathcal{N}$ of measure zero that might depend on $\omega \in \Omega$ (assuring Assumption \ref{as:B-1}).
	Further, the coefficient $a(\omega, x)$ has positive spatially bounded paths, i.e., $0 < a_{-}(\omega) \leq a(\omega, x) \leq a_{+}(\omega) < \infty$ for some constants $a_{-}$ and $a_{+}$ that might depend on $\omega \in \Omega$ (assuring Assumption \ref{as:B-2}).
	The flux function $f (\rho)$ is locally Lipschitz continuous and is either monotone (Assumption \ref{as:B-3alt}) or convex (or concave) with $f(\solution_M (\omega, x)) = 0$ for some $\solution_M$, in which case we have $M_0 = 0$ (Assumption \ref{as:B-3}).
\end{example}

\begin{remark}[Existence of stochastic steady state solutions]
	\label{rem:stochFluxAssumptions}

	To define a meaningful entropy solution, we have to ensure the well-posedness of the stochastic steady state equation
	\begin{equation}
		\label{eq:stochSteadyState}
		\mathfrak{\solFlux}(\omega, x, \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) = \alpha \quad \text{for a.e. }x \in \mathbb{R} \ .
	\end{equation}
	If assumptions \ref{as:B-1}\--\ref{as:B-3} are satisfied, then for any constant $\alpha \in [M_0, \infty)$ (or $\alpha \in (-\infty, M_0]$) and any $\omega \in \Omega$, there exist two steady state solutions $\stochSteadyStateSolution{+}: \Omega \times \mathbb{R} \rightarrow [M_0, \infty)$ and $\stochSteadyStateSolution{-}: \Omega \times \mathbb{R} \rightarrow (- \infty, M_0]$ of \eqref{eq:stochasticConservationLaw} satisfying \eqref{eq:stochSteadyState}.
	For the case \ref{as:B-1}\--\ref{as:B-3alt} the two steady state solutions coincide, i.e., $\stochSteadyStateSolutionArg{+} = \stochSteadyStateSolutionArg{-}$.
\end{remark}

The steady state solutions defined in Remark \ref{rem:stochFluxAssumptions} allow us to introduce the partially adapted Kru\v{z}khov entropies $\partAdaptedEntropies{\alpha}(\omega, x, u) = |u - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)|$, which lead to the notion of pathwise adapted entropy solutions.
\begin{definition}[Pathwise adapted entropy solution]
	Let $\omega \in \Omega$ be fixed.
	A function $u(\omega, \cdot, \cdot) \in L^{\infty}(\mathbb{T} \times \mathbb{R}) \cap C^{0}(\mathbb{T}, L^{1}_{\text{loc}}(\mathbb{R}))$ is an adapted entropy solution of \eqref{eq:stochasticConservationLaw} on $\mathbb{R} \times \mathbb{T}$, provided that for each $\alpha \in [M_0, \infty)$ (or $\alpha \in (-\infty, M_0]$) and the corresponding two steady-state solutions $\entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)$ of \eqref{eq:stochSteadyState}, the inequality
	\begin{equation}
		\label{eq:adaptedEntropyCondition_integralForm}
		\begin{aligned}
			&\int \Big| u(\omega, x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar) \Big| \frac{\partial}{\partial \timeVar} \phi(x, t) \d t \d x \\
			&\hspace*{1.5cm} + \int \sign \Big( u(\omega, x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar) \Big) \Big( \mathfrak{\solFlux}(\omega, x, u(\omega, x, t)) - \alpha \Big) \frac{\partial}{\partial \spaceVar} \phi(x, t) \d t \d x \\
			&\hspace*{3cm} + \int \Big| \solution_{0}(\omega, x) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar) \Big| \phi (x, 0) \d x \geq 0
		\end{aligned}
	\end{equation}
	holds for every test function $\phi \in C^{\infty}_{c}(\mathbb{R} \times \mathbb{T}; \mathbb{R}_{\geq 0})$.
\end{definition}

As a direct consequence of Assumption \ref{as:stochFlux} and the definition of pathwise adapted entropy solutions, we can conclude the following lemma.
\begin{lemma}
	\label{lem:fluxBound}
	Let Assumption \ref{as:B-2} be satisfied and fix an arbitrary $\omega \in \Omega$. Furthermore, suppose $u(\omega, \cdot, \cdot) \in L^{\infty}(\mathbb{T} \times \mathbb{R}) \cap C^{0}(\mathbb{T}, L^{1}_{\text{loc}}(\mathbb{R}))$ is a pathwise adapted entropy solution of \eqref{eq:stochasticConservationLaw}.
	Then the following estimate holds for a.e. $x \in \mathbb{R}$ and a.e. $t \in \mathbb{T}$:
	\begin{equation}
		\label{eq:fluxBound}
		\abs{\mathfrak{\solFlux} (\omega, x, u(\omega, x, t))} \leq \max\limits_{-\norm{u(\omega, \cdot, \cdot)}_{L^{\infty}} \leq \sigma \leq \norm{u(\omega, \cdot, \cdot)}_{L^{\infty}}} g_{+} (\omega, \sigma) =: M(\omega) < \infty \ .
	\end{equation}
\end{lemma}

\begin{proof}
	By hypothesis, we have, for fixed $\omega \in \Omega$, that $u (\omega, \cdot, \cdot) \in L^{\infty}(\mathbb{T} \times \mathbb{R})$.
	Therefore, the maximization is over a compact interval.
	Further, by assumption \ref{as:B-2}, we have the existence of an upper bound $g_{+}$ of $\mathfrak{\solFlux}$, which is continuous in the second argument.
	The maximization is over a compact interval, hence we obtain the existence of the maximum.
	Also, with the same argument, we have $M (\omega) < \infty$, which concludes the proof.
\end{proof}

We are now able to state the following pathwise existence and uniqueness result:
\begin{theorem}
	\label{thm:stochExistenceUniqueness}
	Let the flux function $\mathfrak{\solFlux}$ satisfy Assumption \ref{as:stochFlux}.
	Then, for every $\omega \in \Omega$, there exists a unique pathwise adapted entropy solution $u (\omega, \cdot, \cdot) \in L^{\infty}(\mathbb{T} \times \mathbb{R}) \cap C^{0}(\mathbb{T}, L^{1}_{\text{loc}}(\mathbb{R}))$ of the initial value problem \eqref{eq:stochasticConservationLaw}, which satisfies for a.e. $t \in \mathbb{T}$:
	\begin{equation}
		\label{eq:stochUniqueness}
		\int_{a}^{b} \abs{u(\omega, x, t)} \d{x} \leq \int_{a-M(\omega)t}^{b+M(\omega)t}  \abs{\solution_{0}(\omega, x)} \d{x} \ .
	\end{equation}
	Here, $M (\omega)$ is defined as in \eqref{eq:fluxBound} and $a, b \in \mathbb{R}$ are arbitrary with $a < b$.
\end{theorem}

\begin{proof}
	Let $\omega \in \Omega$ be fixed.
	The existence of a pathwise adapted entropy solution is proved similarly as for the deterministic case, see \cite{Chen2008Hyperbolic}.
	The proof is based on the reduction of measure-valued solutions to entropy solutions in $L^{\infty}$ and the consideration of a mollified version of the problem. \\
	For the uniqueness, which was shown in \cite{Audusse2005Uniqueness}, the use of adapted entropies allows the argumentation via the "doubling of variables" procedure of Kru\v{z}khov \cite{Kruzkov1970First}.
	Then, estimate \eqref{eq:stochUniqueness} is an immediate consequence of the uniqueness result.
\end{proof}

Note, the uniqueness result in Theorem \ref{thm:stochExistenceUniqueness} immediately implies that $u (\omega, \cdot, t) \in L^1(\mathbb{R})$ for fixed $t \in \mathbb{T}$, if $\solution_{0} (\omega, \cdot) \in L^1(\mathbb{R})$.
Furthermore, since we assumed $\solution_{0} \in L^{p}(\Omega; L^{\infty}(\mathbb{R}))$, we immediately obtain the existence of the $p$-th moment of the stochastic entropy solution via Inequality \eqref{eq:stochUniqueness}.

\subsection{Measurability of the solution}

In Theorem \ref{thm:stochExistenceUniqueness}, we have shown the existence and uniqueness of pathwise adapted entropy solutions.
To complete the well-posedness investigation, it remains to show the measurability of the solution map $u: \Omega \rightarrow L^{\infty}(\mathbb{R} \times \timeInterval) \cap C^{0}(\timeInterval; L^{1}_{\text{loc}}(\mathbb{R}))$.
The following assumptions will allow us to prove the measurability of the solution map.
\begin{assumption}[Measurability of the stochastic flux function]
	\label{as:stochFlux_measurability}
	\begin{enumerate}[label=(B-\arabic*)]
		\item[\mylabel{as:B-4}{(B-4)}] For almost every $x \in \mathbb{R}$, the mapping $\omega \mapsto \mathfrak{\solFlux}(\omega, x, \rho)$ is measurable.
		\item[\mylabel{as:B-5}{(B-5)}] For almost every $x \in \mathbb{R}$, the mapping $\omega \mapsto \fluxFunction^{-1}_{\pm}(\omega, x, \rho)$ is measurable, where $\fluxFunction^{-1}_{\pm}(\omega, x, \cdot)$ is the inverse function of $\mathfrak{\solFlux}(\omega, x, \cdot)$.
		\item[\mylabel{as:B-6}{(B-6)}] The spatial paths of $\fluxFunction^{-1}_{\pm}$ are bounded in the sense that, for fixed $\omega \in \Omega, \rho \in \mathbb{R}$, we have $\fluxFunction^{-1}_{\pm}(\omega, \cdot, \rho) \in L^{\infty}(\mathbb{R})$.
	\end{enumerate}
\end{assumption}
\begin{remark}
	\label{rem:measurabilityAssumption}

	The existence of the function $\fluxFunction^{-1}_{\pm}$ is already guaranteed by Assumption \ref{as:stochFlux}, see also Remark \ref{rem:stochFluxAssumptions}.
\end{remark}
\begin{example}[Measurability of multiplicative flux]
	\label{ex:multiplicativeFlux_measurability}

	Considering again the flux defined in Example \ref{ex:multiplicativeFlux}, the Assumption \ref{as:B-4} reduces to the following:
	For almost every $x \in \mathbb{R}$, the mapping $\omega \mapsto a (\omega, x)$ is measurable. \\
	For Assumption \ref{as:B-5}, we have
	\begin{equation}
		\fluxFunction^{-1}_{\pm}(\omega, x, \rho) = \frac{1}{a(\omega, x)} f^{-1}(\rho) .
	\end{equation}
	Since $f^{-1}(\rho)$ is continuous, it is especially measurable and thus, $\fluxFunction^{-1}_{\pm}(\omega, x, \rho)$ is measurable as the composition of two measurable functions.
	Consequently, for multiplicative flux functions satisfying Assumptions \ref{as:B-1}\--\ref{as:B-2} and \ref{as:B-3} or \ref{as:B-3alt}, together with Assumption \ref{as:B-4}, Assumption \ref{as:B-5} is automatically satisfied.
	Finally, since we have $0 < a_{-}(\omega) \leq a(\omega, x) \leq a_{+}(\omega) < \infty$, the spatial paths of $\fluxFunction^{-1}_{\pm}$ are bounded.
\end{example}

To show the measurability of the solution map $u: \Omega \rightarrow L^{\infty}(\mathbb{R} \times \timeInterval) \cap C^{0}(\timeInterval; L^{1}_{\text{loc}}(\mathbb{R}))$, we will need the following \textit{adapted entropy functional}.
For the sake of notational simplicity, we define the function space $\mathcal{X} \coloneqq L^{\infty}(\mathbb{R} \times \timeInterval) \cap C^{0}(\timeInterval; L^{1}_{\text{loc}}(\mathbb{R}))$.
Note that $\mathcal{X}$ is separable, since $L^{1}_{\text{loc}}(\mathbb{R})$ is separable (due to $L^{1}(\mathbb{R})$ being separable and a dense subspace of $L^{1}_{\text{loc}}(\mathbb{R})$) and $\mathbb{T}$ is compact.
\begin{definition}[Adapted entropy functional]
	\label{def:entropyFunctional}

	Let Assumption \ref{as:stochFlux} be satisfied.
	Further, let $\alpha \in [M_0, \infty)$ (or $(-\infty, M_0]$) be fixed with corresponding steady state solutions $\stochSteadyStateSolution{\pm}$ defined by Equation \eqref{eq:stochSteadyState}.
	Then, the adapted entropy functional $J^{\alpha}: \Omega \times \mathcal{X} \rightarrow \mathbb{R}$ is defined as
	\begin{subequations}
		\label{eq:entropyFunctional}
		\begin{align}
			(\omega, \nu) &\mapsto \int \abs{\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} \frac{\partial}{\partial \timeVar} \phi(x, t) \d t \d x \label{eq:entropyIntegral-1}\\
			&\qquad+ \int \sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \frac{\partial}{\partial \spaceVar} \phi(x, t) \d t \d x \label{eq:entropyIntegral-2} \\
			&\qquad\qquad + \int \abs{\solution_{0}(\omega, x) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} \phi (x, 0) \d x \label{eq:entropyIntegral-3} ,
		\end{align}
	\end{subequations}
	for some test function $\phi \in C^{\infty}_{c}(\mathbb{R} \times \mathbb{T}; \mathbb{R}_{\geq 0})$.
\end{definition}
\begin{remark}[Adapted entropy condition via adapted entropy functional]
	Note, with the adapted entropy functional, the adapted entropy condition \eqref{eq:adaptedEntropyCondition_integralForm} can be formulated as: \\
	A function $u (\omega, \cdot, \cdot)$ is an adapted entropy solution, provided that for each $\alpha \in [M_0, \infty)$ (or $(-\infty, M_0]$) and the corresponding two steady state solutions $\entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)$, it holds that $J^{\alpha}(\omega, u) \geq 0$, for every test function $\phi \in C^{\infty}_{c}(\mathbb{R} \times \mathbb{T}; \mathbb{R}_{\geq 0})$.
\end{remark}

Before we are able to prove the measurability of the solution map, we need to establish a few properties of the steady state equations $\stochSteadyStateSolutionArg{\pm}$ and the adapted entropy functional $J^{\alpha}$.
\begin{corollary}[Measurability of steady state solutions]
	\label{cor:measurability_steadyState}
	Let Assumptions \ref{as:stochFlux} and \ref{as:stochFlux_measurability} be satisfied. Then, for any constant $\alpha \in [M_0, \infty)$ (or $\alpha \in (-\infty, M_0]$), the stochastic steady state solutions $\stochSteadyStateSolution{\pm}$ defined by Equation \eqref{eq:stochSteadyState} are measurable in the sense that, for almost every $x \in \mathbb{R}$, the mapping $\omega \mapsto \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)$ is measurable.
\end{corollary}

\begin{proof}
	By assumptions \ref{as:B-1}\--\ref{as:B-3}, we have the existence of the steady state solutions, which are measurable due to Assumption \ref{as:B-5}.
	Here, we mean measurability in the sense that, for each $x \in \mathbb{R}$, the mapping $\omega \mapsto \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)$ is measurable.
\end{proof}
\begin{lemma}[Steady state solutions depend continuously on $\alpha$]
	\label{lem:steadyStateContinuity}
	Let Assumption \ref{as:stochFlux} be satisfied.
	Then, the steady state solutions defined via Equation \eqref{eq:stochSteadyState} depend continuously on $\alpha \in [M_0, \infty)$ (or $\alpha \in (-\infty, M_0]$).
\end{lemma}

\begin{proof}
	Let $\alpha, \beta \in [M_0, \infty)$ (or $\alpha, \beta \in (-\infty, M_0]$) with $\abs{\alpha - \beta} < \delta$.
	Consider the two steady state equations
	\begin{subequations}
		\begin{align}
			\mathfrak{\solFlux}(\omega, x, m^{\pm}_{\alpha}) &= \alpha \quad \text{for a.e. } x \in \mathbb{R} \label{eq:steadyStateAlpha} , \\
			\mathfrak{\solFlux}(\omega, x, m^{\pm}_{\beta}) &= \beta \quad \text{for a.e. } x \in \mathbb{R} \label{eq:steadyStateBeta} .
		\end{align}
	\end{subequations}
	By assumption, we know that there exist two steady state solutions $m^{\pm}_{\alpha}$ and $m^{\pm}_{\beta}$ for each \eqref{eq:steadyStateAlpha} and \eqref{eq:steadyStateBeta}, respectively.
	Thus, we may rewrite these equations as
	\begin{subequations}
		\begin{align}
			m^{\pm}_{\alpha} &= \fluxFunction^{-1}_{\pm}(\omega, x, \alpha) \quad \text{for a.e. } x \in \mathbb{R} \label{eq:steadyStateInvAlpha} , \\
			m^{\pm}_{\beta} &= \fluxFunction^{-1}_{\pm}(\omega, x, \beta) \quad \text{for a.e. } x \in \mathbb{R} \label{eq:steadyStateInvBeta} .
		\end{align}
	\end{subequations}
	Here $\fluxFunction^{-1}_{\pm}(\omega, x, \cdot)$ denotes the inverse of $\mathfrak{\solFlux}(\omega, x, \cdot)$.
	Further, since $\mathfrak{\solFlux}(\omega, x, \cdot)$ is a continuous map from $\mathbb{R}$ to $[M_0, \infty)$ (or $(-\infty, M_0]$), we get that $\fluxFunction^{-1}_{\pm}(\omega, x, \cdot)$ is a continuous function from either $[M_0, \infty)$ or $(-\infty, M_0]$ to $\mathbb{R}$.
	Now, consider
	\begin{equation}
		\begin{aligned}
			\abs{m^{\pm}_{\alpha} - m^{\pm}_{\beta}} = \abs{\fluxFunction^{-1}_{\pm}(\omega, x, \alpha) - \fluxFunction^{-1}_{\pm}(\omega, x, \beta)}
			\leq C_{\omega} \abs{\alpha - \beta} ,
		\end{aligned}
	\end{equation}
	where we used the continuity of $\fluxFunction^{-1}_{\pm}(\omega, x, \cdot)$.
	Note, by Assumption \ref{as:B-6} the constant $C_{\omega}$ does not depend on $x \in \mathbb{R}$.
	Consequently, we have shown that the steady state solutions \eqref{eq:stochSteadyState} depend continuously on $\alpha$.
\end{proof}
\begin{proposition}[Adapted entropy functional depends continuously on $\alpha$]
	Let the flux function satisfy Assumptions \ref{as:stochFlux} and \ref{as:stochFlux_measurability}.
	Let $\alpha \in [M_0, \infty)$ (or $(-\infty, M_0]$) be fixed and let $\stochSteadyStateSolution{\pm}$ denote the steady state solutions defined by Equation \eqref{eq:stochSteadyState}.
	Then the adapted entropy functional $J^{\alpha}$ defined in Equation \eqref{eq:entropyFunctional} depends continuously on $\alpha$.
\end{proposition}

\begin{proof}
	Let $\alpha, \beta \in [M_0, \infty)$ (or $(-\infty, M_0]$) with $\abs{\alpha - \beta} < \delta$.
	Due to Lemma \ref{lem:steadyStateContinuity}, it immediately follows that the terms \eqref{eq:entropyIntegral-1} and \eqref{eq:entropyIntegral-3} are continuous in $\alpha$.
	To show the continuity of \eqref{eq:entropyIntegral-2}, we define
	\begin{equation}
		\label{eq:G_defintion}
		G^{\alpha}(\omega, \nu) \coloneqq \int \sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \frac{\partial}{\partial \spaceVar} \phi(x, t) \d t \d x ,
	\end{equation}
	and consider
	\begin{equation}
		\label{eq:continuityDef_alpha}
		\begin{aligned}
			\abs{ G^{\alpha} (\omega, \nu(x, t)) - G^{\beta}(\omega, \nu(x, t))} &= \Big| \int \Big( \sign (\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&\hspace*{2cm} - \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta ) \Big) \frac{\partial}{\partial \spaceVar} \phi \d t \d x \Big| .
		\end{aligned}
	\end{equation}

	Define the following three sets:
	\begin{subequations}
		\begin{align}
			D^{=} &= \set{x \in \mathbb{R} : \sign(\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) = \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x))} , \\
			D^{\pm} &= \set{x \in \mathbb{R} : \sign(\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) = - \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x))} , \\
			D^{m} &= \set{x \in \mathbb{R} : \nu(x, t) = m^{\pm}_{\alpha}(\omega, x) \text{ or } \nu(x, t) = m^{\pm}_{\beta}(\omega, x)} .
		\end{align}
	\end{subequations}

	Note, these sets are disjoint, whenever we have $\alpha \neq \beta$, since we have
	\begin{equation}
		D^{=} \cap D^{\pm} \cap D^{m} = \set{x \in \mathbb{R}: m^{\pm}_{\alpha}(\omega, x) = \nu(x, t) = m^{\pm}_{\beta}(\omega, x)} .
	\end{equation}
	Hence, we can rewrite Equation \eqref{eq:continuityDef_alpha} as
	\begin{subequations}
		\label{eq:integralSplit_alpha}
		\begin{align}
			\Big| G^{\alpha} (\omega, \nu(x, t)) &- G^{\beta}(\omega, \nu(x, t)) \Big|  \\
			\begin{split}
				&\leq \Big|\int_{D^{=}} \Big( \sign (\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
				&\hspace*{3cm} - \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta ) \Big) \frac{\partial}{\partial \spaceVar} \phi \d t \d x \Big| \label{eq:integral_D=_alpha}
			\end{split}\\
			\begin{split}
				&\qquad + \Big| \int_{D^{\pm}} \Big( \sign (\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
				&\hspace*{3cm} - \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta ) \Big) \frac{\partial}{\partial \spaceVar} \phi \d t \d x \Big| \label{eq:integral_Dpm_alpha}
			\end{split}\\
			\begin{split}
				&\qquad + \Big| \int_{D^{m}} \Big( \sign (\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
				&\hspace*{3cm} - \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta ) \Big) \frac{\partial}{\partial \spaceVar} \phi \d t \d x \Big| \label{eq:integral_Dm_alpha} .
			\end{split}
		\end{align}
	\end{subequations}

	We will now estimate the terms in \eqref{eq:integralSplit_alpha} separately, starting with \eqref{eq:integral_D=_alpha}:
	\begin{equation}
		\begin{aligned}
			\Big| \int_{D^{=}} \Big( \sign (\nu(x, t) &- m^{\pm}_{\alpha}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&- \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta ) \Big) \frac{\partial}{\partial \spaceVar} \phi \d t \d x \Big| \\
			&\hspace*{2cm}\leq \int_{D^{=}} \Big| \sign (\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&\hspace*{4cm} - \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta ) \frac{\partial}{\partial \spaceVar} \phi \Big| \d t \d x \\
			&\hspace*{2cm}\leq \int_{D^{=}} \abs{ \Big( \mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha - \mathfrak{\solFlux}(\omega, x, \nu(x, t)) + \beta \Big) \frac{\partial}{\partial \spaceVar} \phi} \d t \d x \\
			&\hspace*{2cm}\leq C_{\phi} \abs{\alpha - \beta} .
		\end{aligned}
	\end{equation}
	Here, we first used the definition of $D^{=}$ and afterwards the boundedness of $\frac{\partial}{\partial \spaceVar} \phi$.

	Now, consider the integral \eqref{eq:integral_Dpm_alpha}:
	\begin{equation}
		\begin{aligned}
			\Big| \int_{D^{\pm}} \Big( \sign (&\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&- \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta ) \Big) \frac{\partial}{\partial \spaceVar} \phi \d t \d x \Big| \\
			&\hspace*{1.5cm}\leq \int_{D^{\pm}} \Big| \sign (\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&\hspace*{3cm} - \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta ) \frac{\partial}{\partial \spaceVar} \phi \Big| \d t \d x \\
			&\hspace*{1.5cm}\leq \int_{D^{\pm}} \abs{ (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \frac{\partial}{\partial \spaceVar} \phi } \d t \d x + \int_{D^{\pm}} \abs{ (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta) \frac{\partial}{\partial \spaceVar} \phi } \d t \d x \\
			&\hspace*{1.5cm}\leq \int_{D^{\pm}} L_{I} \abs{ \nu(x, t) - m^{\pm}_{\alpha}(\omega, x)} \abs{ \frac{\partial}{\partial \spaceVar} \phi } \d t \d x + \int_{D^{\pm}} L_{I} \abs{ \nu(x, t) - m^{\pm}_{\beta}(\omega, x)} \abs{ \frac{\partial}{\partial \spaceVar} \phi } \d t \d x .
		\end{aligned}
	\end{equation}
	In the second step, we used the definition of $D^{\pm}$ together with the triangle inequality and estimated $\sign(\nu - m^{\pm}_{\alpha}) \leq 1$.
	Afterwards, we inserted the definition of $\alpha$ and $\beta$ via Equations \eqref{eq:steadyStateAlpha} and \eqref{eq:steadyStateBeta}, respectively, and used the local Lipschitz continuity of $\mathfrak{\solFlux}$.
	Using now the fact that for $\xi_1, \xi_2 \in \mathbb{R}$, we have the identity $\abs{\xi_1 - \xi_2} = \sign(\xi_1 - \xi_2)(\xi_1 - \xi_2)$, we obtain
	\begin{equation}
		\begin{aligned}
			\Big| \int_{D^{\pm}} \Big( \sign (&\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&- \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta ) \Big) \frac{\partial}{\partial \spaceVar} \phi \d t \d x \Big| \\
			&\hspace*{0.5cm}\leq \int_{D^{\pm}} L_{I} \Big( \sign(\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) (\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) \\
			&\hspace*{3cm} + \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x)) (\nu(x, t) - m^{\pm}_{\beta}(\omega, x)) \Big) \abs{ \frac{\partial}{\partial \spaceVar} \phi } \d t \d x \\
			&\hspace*{0.5cm}= \int_{D^{\pm}} L_{I} \sign(\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) \Big( \nu(x, t) - m^{\pm}_{\alpha}(\omega, x) - \nu(x, t) + m^{\pm}_{\beta}(\omega, x) \Big) \abs{ \frac{\partial}{\partial \spaceVar} \phi } \d t \d x \\
			&\hspace*{0.5cm}\leq \int_{D^{\pm}} L_{I} \abs{ m^{\pm}_{\alpha}(\omega, x) - m^{\pm}_{\beta}(\omega, x)} \abs{\frac{\partial}{\partial \spaceVar} \phi} \d t \d x \\
			&\hspace*{0.5cm}\leq L_{I} C_{\phi} \abs{ m^{\pm}_{\alpha}(\omega, x) - m^{\pm}_{\beta}(\omega, x)} \\
			&\hspace*{0.5cm}\leq L_{I} C_{\phi} C_{\omega} \abs{\alpha - \beta}.
		\end{aligned}
	\end{equation}
	Here, we used in the last step the continuity of $m^{\pm}_{\alpha}$ in $\alpha$ from Lemma \ref{lem:steadyStateContinuity}.
	It remains to estimate \eqref{eq:integral_Dm_alpha}.
	Therefore, we define the two sets
	\begin{subequations}
		\begin{align}
			D^{m}_{\alpha} &= \set{x \in \mathbb{R}: \nu(x, t) = m^{\pm}_{\alpha}(\omega, x)} , \\
			D^{m}_{\beta} &= \set{x \in \mathbb{R}: \nu(x, t) = m^{\pm}_{\beta}(\omega, x)} .
		\end{align}
	\end{subequations}
	Note that these sets are disjoint, whenever $\alpha \neq \beta$.
	With the definition of $D^{m}_{\alpha}$ and $D^{m}_{\beta}$ and the triangle inequality, we obtain for \eqref{eq:integral_Dm_alpha}
	\begin{equation}
		\begin{aligned}
			\Big| \int_{D^{m}} \Big( \sign (\nu(x, t) &- m^{\pm}_{\alpha}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&- \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta ) \Big) \frac{\partial}{\partial \spaceVar} \phi \d t \d x \Big| \\
			&\hspace*{1.5cm}\leq \Big| \int_{D^{m}_{\alpha}} - \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta ) \frac{\partial}{\partial \spaceVar} \phi \d t \d x \Big| \\
			&\hspace*{3cm} + \Big| \int_{D^{m}_{\beta}} \sign (\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \frac{\partial}{\partial \spaceVar} \phi \d t \d x \Big| \\
			&\hspace*{1.5cm}\leq \int_{D^{m}_{\alpha}} \abs{\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta} \abs{\frac{\partial}{\partial \spaceVar} \phi} \d t \d x + \int_{D^{m}_{\beta}} \abs{\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha} \abs{\frac{\partial}{\partial \spaceVar} \phi} \d t \d x . \\
		\end{aligned}
	\end{equation}
	Now, using the steady state equations \eqref{eq:steadyStateAlpha} and \eqref{eq:steadyStateBeta} and the local Lipschitz continuity of $\mathfrak{\solFlux}$, we get
	\begin{equation}
		\begin{aligned}
			\Big| \int_{D^{m}} \Big( &\sign (\nu(x, t) - m^{\pm}_{\alpha}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&- \sign(\nu(x, t) - m^{\pm}_{\beta}(\omega, x)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \beta ) \Big) \frac{\partial}{\partial \spaceVar} \phi \d t \d x \Big| \\
			&\hspace*{1.5cm}= \int_{D^{m}_{\alpha}} \abs{\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \mathfrak{\solFlux}(\omega, x, m^{\pm}_{\beta}(\omega, x))} \abs{\frac{\partial}{\partial \spaceVar} \phi} \d t \d x \\
			&\hspace*{3cm}+ \int_{D^{m}_{\beta}} \abs{\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \mathfrak{\solFlux}(\omega, x, m^{\pm}_{\alpha}(\omega, x))} \abs{\frac{\partial}{\partial \spaceVar} \phi} \d t \d x \\
			&\hspace*{1.5cm}= \int_{D^{m}_{\alpha}} \abs{\mathfrak{\solFlux}(\omega, x, m^{\pm}_{\alpha}(\omega, x)) - \mathfrak{\solFlux}(\omega, x, m^{\pm}_{\beta}(\omega, x))} \abs{\frac{\partial}{\partial \spaceVar} \phi} \d t \d x \\
			&\hspace*{3cm}+ \int_{D^{m}_{\beta}} \abs{\mathfrak{\solFlux}(\omega, x, m^{\pm}_{\beta}(\omega, x)) - \mathfrak{\solFlux}(\omega, x, m^{\pm}_{\alpha}(\omega, x))} \abs{\frac{\partial}{\partial \spaceVar} \phi} \d t \d x \\
			&\hspace*{1.5cm}\leq \int_{D^{m}_{\alpha}} L_{I} \abs{m^{\pm}_{\alpha}(\omega, x) - m^{\pm}_{\beta}(\omega, x)} \abs{\frac{\partial}{\partial \spaceVar} \phi} \d t \d x + \int_{D^{m}_{\beta}} L_{I} \abs{m^{\pm}_{\beta}(\omega, x) - m^{\pm}_{\alpha}(\omega, x)} \abs{\frac{\partial}{\partial \spaceVar} \phi} \d t \d x \\
			&\hspace*{1.5cm}\leq 2 L_{I} C_{\phi} C_{\omega} \abs{\alpha - \beta} .
		\end{aligned}
	\end{equation}
	Here, we used Lemma \ref{lem:steadyStateContinuity} in the last step.
	Consequently, we get the following estimation for \eqref{eq:integralSplit_alpha}:
	\begin{equation}
		\begin{aligned}
			\Big| G^{\alpha} (\omega, \nu) - G^{\beta}(\omega, \nu) \Big| &\leq C_{\phi} \abs{\alpha - \beta} + L_{I} C_{\phi} C_{\omega} \abs{\alpha - \beta} + 2 L_{I} C_{\phi} C_{\omega} \abs{\alpha - \beta} \\
			&\leq \Big( C_{\phi} + 3 L_{I} C_{\phi} C_{\omega} \Big) \abs{\alpha - \beta}
		\end{aligned}
	\end{equation}
	and thus, the functional $J^{\alpha}$ depends continuously on $\alpha$.
\end{proof}

As a last property of the adapted entropy functional, we will show that $J^{\alpha}$ is a Carathéodory functional.
\begin{proposition}[Adapted entropy functional is Carathéodory]
	\label{prop:EntropyFunctional_Caratheodory}
	Let the flux function $\mathfrak{\solFlux}$ satisfy Assumptions \ref{as:stochFlux} and \ref{as:stochFlux_measurability} and let $\solution_{0} \in L^{p}(\Omega, L^{\infty}(\mathbb{R}))$.
	Further, let $\alpha \in [M_0, \infty)$ (or $(-\infty, M_0]$) be fixed and let $\stochSteadyStateSolution{\pm}$ denote the steady state solutions defined by Equation \eqref{eq:stochSteadyState}.
	Then the adapted entropy functional $J^{\alpha}$, defined in \eqref{eq:entropyFunctional} is a Carathéodory map, i.e., $J^{\alpha}(\omega, \nu)$ is measurable in $\omega \in \Omega$ and continuous in $\nu \in \mathcal{X}$.
\end{proposition}

\begin{proof}
	We divide the step into two steps: we first show the continuity with respect to $\nu \in \mathcal{X}$ and afterwards the measurability in $\omega \in \Omega$.

	\textit{Step 1:}
	We start by showing the continuity of $J^{\alpha}$ w.r.t. $\nu \in \mathcal{X}$.
	The continuity of the first integral \eqref{eq:entropyIntegral-1} and the third integral \eqref{eq:entropyIntegral-3} is obvious.
	For the second integral \eqref{eq:entropyIntegral-2}, we explicitly show the continuity: \\
	Let $\omega \in \Omega$ and $\alpha \in [M_0, \infty)$ (or $(-\infty, M_0]$) be fixed and let $\nu, \tilde{\nu} \in \mathcal{X}$ with $\norm{\nu - \tilde{\nu}}_{\mathcal{X}} < \delta$, for some $\delta > 0$.
	Consider now
	\begin{equation}
		\begin{aligned}
			\abs{G^{\alpha}(\omega, \nu) - G^{\alpha}(\omega, \tilde{\nu})} &= \Big| \int \Big(\sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&\qquad \qquad \qquad - \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha) \Big) \frac{\partial}{\partial \spaceVar} \phi(x, t) \d t \d x \Big| ,
		\end{aligned}
	\end{equation}
	where $G^{\alpha}(\omega, \nu)$ is defined as in Equation \eqref{eq:G_defintion}.
	Let us define the following three sets:
	\begin{equation}
		\begin{aligned}
			D^{=} &\coloneqq \set{x \in \mathbb{R} | \sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) = \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar))} , \\
			D^{m} &\coloneqq \set{x \in \mathbb{R} | \nu(x, t) = \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar) \text{ or } \tilde{\nu}(x, t) = \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} , \\
			D^{\pm} &\coloneqq \set{x \in \mathbb{R} | \sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) = - \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar))} .
		\end{aligned}
	\end{equation}

	Note that in general we have $D^{=} \cap D^{m} \cap D^{\pm} \neq \emptyset$, but we have $\mathbb{R} = D^{=} \cup D^{m} \cup D^{\pm}$.
	Thus, we can write
	\begin{subequations}
		\label{eq:continuityEstimate}
		\begin{align}
			\abs{G^{\alpha}(\omega, \nu) - G^{\alpha}(\omega, \tilde{\nu})} &= \Big| \int \Big(\sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&\qquad \qquad \qquad - \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha) \Big) \frac{\partial}{\partial \spaceVar} \phi(x, t) \d t \d x \Big| \\
			\begin{split}
				&\leq \int_{D^{=} \times \mathbb{T}} \Big| \sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \label{eq:integral_D=} \\
				&\qquad \qquad \qquad - \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x, t)} \d t \d x
			\end{split}\\
			\begin{split}
				&\qquad + \int_{D^{m} \times \mathbb{T}} \Big| \sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \label{eq:integral_Dm} \\
				&\qquad \qquad \qquad - \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x, t)} \d t \d x
			\end{split}\\
			\begin{split}
				&\qquad + \int_{D^{\pm} \times \mathbb{T}} \Big| \sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \label{eq:integral_Dpm} \\
				&\qquad \qquad \qquad - \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x, t)} \d t \d x.
			\end{split}
		\end{align}
	\end{subequations}

	When considering the integral \eqref{eq:integral_D=}, we obtain
	\begin{equation}
		\begin{aligned}
			\int\limits_{D^{=} \times \mathbb{T}} \Big| \Big( \sign &(\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			& - \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha) \Big) \frac{\partial}{\partial \spaceVar} \phi(x, t) \Big| \d t \d x \\[5pt]
			&\hspace*{3cm} = \int\limits_{D^{=} \times \mathbb{T}} \Big| \Big( (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) - (\mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha) \Big) \frac{\partial}{\partial \spaceVar} \phi(x, t) \Big| \d t \d x \\
			&\hspace*{3cm} = \int\limits_{D^{=} \times \mathbb{T}} \Big| \mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x, t)} \d t \d x \\
			&\hspace*{3cm} \leq \int\limits_{D^{=} \times \mathbb{T}} L_I \abs{\nu(x, t) - \tilde{\nu}(x, t)} \abs{\frac{\partial}{\partial \spaceVar} \phi(x, t)} \d t \d x \\
			&\hspace*{3cm} \leq C_{\phi} L_I \norm{\nu - \tilde{\nu}}_{L^{\infty}(D^{=} \times \mathbb{T})} .
		\end{aligned}
	\end{equation}
	Here, we first used the definition of $D^{=}$ to estimate the sign functions.
	Afterwards, we used the locally Lipschitz continuity of $\mathfrak{\solFlux}$ in the third argument by Assumption \ref{as:B-3} or \ref{as:B-3alt}.
	Note that the Lipschitz constant depends on the interval $I \subset \mathbb{R}$, which contains $\nu(x,t), \tilde{\nu}(x,t)$ for all $(x,t) \in \mathbb{R} \times \mathbb{T}$.
	This interval exists, since by hypothesis we have $\nu, \tilde{\nu} \in \mathcal{X}$, which implies $\nu, \tilde{\nu} \in L^{\infty}(\mathbb{R} \times \mathbb{T})$.

	Now, before we consider the integral \eqref{eq:integral_Dm}, we define
	\begin{equation}
		\label{eq:Dm_splitting}
		\begin{aligned}
			D^{m}_{\nu} &\coloneqq \set{x \in D^{m} | \nu(x, t) = \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} \\
			D^{m}_{\tilde{\nu}} &\coloneqq \set{x \in D^{m} | \tilde{\nu}(x, t) = \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} .
		\end{aligned}
	\end{equation}
	With this, we can rewrite the integral \eqref{eq:integral_Dm} as
	\begin{equation}
		\begin{aligned}
			\int\limits_{D^{m} \times \mathbb{T}} \Big| \sign (&\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&- \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x, t)} \d t \d x \\
			&\hspace*{1cm} = \int\limits_{D^{m}_{\nu} \times \mathbb{T}} \Big| \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x, t)} \d t \d x \\
			&\hspace*{3cm} + \int\limits_{D^{m}_{\tilde{\nu}} \times \mathbb{T}} \Big| \sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x, t)} \d t \d x \\
			&\hspace*{1cm} \leq \int\limits_{D^{m}_{\nu} \times \mathbb{T}} \Big|  \mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha \Big| \abs{ \frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x + \int\limits_{D^{m}_{\tilde{\nu}} \times \mathbb{T}} \Big| \mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x, t)} \d t \d x .
		\end{aligned}
	\end{equation}
	Now, note that Equation \eqref{eq:Dm_splitting} gives by definition
	\begin{equation}
		\begin{aligned}
			\alpha &= \mathfrak{\solFlux}(\omega, x, \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) = \mathfrak{\solFlux}(\omega, x, \nu(x, t)) \qquad \text{for } x \in D^{m}_{\nu} , \\
			\alpha &= \mathfrak{\solFlux}(\omega, x, \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) = \mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) \qquad \text{for } x \in D^{m}_{\tilde{\nu}} .
		\end{aligned}
	\end{equation}
	Consequently, we obtain
	\begin{equation}
		\begin{aligned}
			\int\limits_{D^{m} \times \mathbb{T}} \Big| \sign (&\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&- \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x, t)} \d t \d x \\
			&\hspace*{1cm} \leq \int\limits_{D^{m}_{\nu} \times \mathbb{T}} \Big|  \mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha \Big| \abs{ \frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x + \int\limits_{D^{m}_{\tilde{\nu}} \times \mathbb{T}} \Big| \mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x, t)} \d t \d x \\
			&\hspace*{1cm} = \int_{D^{m}_{\nu} \times \mathbb{T}} \Big|  \mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \mathfrak{\solFlux}(\omega, x, \nu(x, t)) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{3cm} + \int_{D^{m}_{\tilde{\nu}} \times \mathbb{T}} \Big| \mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{1cm} \leq C_{\phi} L_I \norm{\nu - \tilde{\nu}}_{L^{\infty}(D^{m}_{\nu} \times \mathbb{T})}  + C_{\phi} L_I \norm{\nu - \tilde{\nu}}_{L^{\infty}(D^{m}_{\tilde{\nu}} \times \mathbb{T})} \\
			&\hspace*{1cm} \leq 2 C_{\phi} L_I \norm{\nu - \tilde{\nu}}_{L^{\infty}(D^{m} \times \mathbb{T})} .
		\end{aligned}
	\end{equation}
	It remains to consider the integral \eqref{eq:integral_Dpm}:
	\begin{equation}
		\begin{aligned}
			\int_{D^{\pm} \times \mathbb{T}} \Big| \sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(&\omega, x, \nu(x, t)) - \alpha) \\
			&- \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x.
		\end{aligned}
	\end{equation}

	We start with the estimation
	\begin{equation}
		\begin{aligned}
			\int\limits_{D^{\pm} \times \mathbb{T}} \Big| \sign (&\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&- \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{2cm} \leq \int\limits_{D^{\pm} \times \mathbb{T}} \Big| \mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{4cm} + \int\limits_{D^{\pm} \times \mathbb{T}} \Big| \mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{2cm} = \int\limits_{D^{\pm} \times \mathbb{T}} \Big| \mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \mathfrak{\solFlux}(\omega, x, \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{4cm} + \int\limits_{D^{\pm} \times \mathbb{T}} \Big| \mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \mathfrak{\solFlux}(\omega, x, \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{2cm} \leq \int\limits_{D^{\pm} \times \mathbb{T}} L_I \abs{\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{4cm} + \int\limits_{D^{\pm} \times \mathbb{T}} L_I \abs{\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{2cm} \leq \int\limits_{D^{\pm} \times \mathbb{T}} L_I \Big( \abs{\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} + \abs{\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} \Big) \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x . \\
		\end{aligned}
	\end{equation}
	Here, we first used the steady state equation \eqref{eq:stochSteadyState} and then the local Lipschitz continuity of $\mathfrak{\solFlux}$ in the third argument.
	Note that the Lipschitz constant is dependent on the interval $I$ containing $\nu(x,t), \tilde{\nu}(x,t)$ and $\entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)$ for all $(x,t) \in \mathbb{R} \times \mathbb{T}$.
	The interval exists, since $\nu, \tilde{\nu} \in L^{\infty}(\mathbb{R} \times \mathbb{T})$ by hypothesis and, for each $\omega, \alpha$, we have $\entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar) \in L^{\infty}(\mathbb{R} \times \mathbb{T})$ due to Assumption \ref{as:B-2}.
	Note, even though we do not explicitly know that the above integrals are well defined, their existence is ensured via dominated convergence with the subsequent estimation.

	By the definition of $D^{\pm}$, we have
	\begin{equation}
		\label{eq:D_plusminus}
		\sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) = - \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) \qquad \text{for } x \in D^{\pm} .
	\end{equation}

	For $\xi_1, \xi_2 \in \mathbb{R}$ it holds that $\abs{\xi_1 - \xi_2} = \sign(\xi_1 - \xi_2)(\xi_1 - \xi_2)$.
	Together with \eqref{eq:D_plusminus} this yields
	\begin{equation}
		\begin{aligned}
			\int\limits_{D^{\pm} \times \mathbb{T}} \Big| \sign (&\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \\
			&- \sign (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \tilde{\nu}(x, t)) - \alpha) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{2cm} \leq \int\limits_{D^{\pm} \times \mathbb{T}} L_I \Big( \abs{\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} + \abs{\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} \Big) \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{2cm} = \int\limits_{D^{\pm} \times \mathbb{T}} L_I \Big( \sign(\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar))(\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) \\
			&\hspace*{4cm} + \sign(\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar))(\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) \Big) \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{2cm} = \int\limits_{D^{\pm} \times \mathbb{T}} L_I \Big( \sign(\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar))(\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) \\
			&\hspace*{4cm}  - \sign(\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar))(\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) \Big) \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{2cm} = L_I \int\limits_{D^{\pm} \times \mathbb{T}} \sign(\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) \\
			&\hspace*{4cm} \cdot \Big( (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) - (\tilde{\nu}(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) \Big) \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{2cm} = L_I \int\limits_{D^{\pm} \times \mathbb{T}} \sign(\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) \Big( \nu(x, t) - \tilde{\nu}(x, t) \Big) \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{2cm} \leq L_I \int\limits_{D^{\pm} \times \mathbb{T}} \Big| \nu(x, t) - \tilde{\nu}(x, t) \Big| \abs{\frac{\partial}{\partial \spaceVar} \phi(x,t)} \d t \d x \\
			&\hspace*{2cm} \leq C_{\phi} L_I \norm{\nu - \tilde{\nu}}_{L^{\infty}(D^{\pm} \times \mathbb{T})}.
		\end{aligned}
	\end{equation}

	Consequently, the continuity estimation \eqref{eq:continuityEstimate} reduces to
	\begin{equation}
		\begin{aligned}
			\abs{G^{\alpha}(\omega, \nu) - G^{\alpha}(\omega, \tilde{\nu})} &\leq C_{\phi} L_I \norm{\nu - \tilde{\nu}}_{L^{\infty}(D^{=} \times \mathbb{T})} + 2 C_{\phi} L_I \norm{\nu - \tilde{\nu}}_{L^{\infty}(D^{m} \times \mathbb{T})} + C_{\phi} L_I \norm{\nu - \tilde{\nu}}_{L^{\infty}(D^{\pm} \times \mathbb{T})} \\
			&\leq 4 C_{\phi} L_I \norm{\nu - \tilde{\nu}}_{L^{\infty}(\mathbb{R} \times \mathbb{T})} ,
		\end{aligned}
	\end{equation}
	and thus,
	\begin{equation}
		G^{\alpha}(\omega, \nu) \coloneqq \int \sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \frac{\partial}{\partial \spaceVar} \phi \d t \d x
	\end{equation}
	is continuous in $\nu \in \mathcal{X}$.

	\textit{Step 2:}
	It remains to show the measurability of $J^{\alpha}$ with resprect to $\omega \in \Omega$.
	Therefore, we first establish that all integrands in Equation \eqref{eq:entropyFunctional} are measurable: \\
	Let $\alpha \in [M_0, \infty)$ (or $(-\infty, M_0]$) be arbitrary but fixed.
	Note that by Corollary \ref{cor:measurability_steadyState}, the steady state solutions $\entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)$ are measurable as a mapping $\omega \mapsto \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)$.
	Hence, the integrand of \eqref{eq:entropyIntegral-1}, i.e., $\abs{\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} \frac{\partial}{\partial \timeVar} \phi(x, t)$ is measurable.
	Since $\entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)$ is measurable, the measurability of the integrand of \eqref{eq:entropyIntegral-2} is an immediate consequence of Assumption \ref{as:B-4}.
	Finally, the integrand of \eqref{eq:entropyIntegral-3} is measurable, since the steady state solutions $\entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)$ are measurable and $\solution_{0} \in L^{p}(\Omega, L^{\infty}(\mathbb{R}))$ is measurable by assumption.

	Now, since $\phi \in C^{\infty}_{c}(\mathbb{R} \times \mathbb{T})$ has compact support we can write $J^{\alpha}$ as
	\begin{subequations}
		\label{eq:functional_compactSet}
		\begin{align}
			(\omega, \nu) &\mapsto \int\limits_{\supp{\phi}} \abs{\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} \frac{\partial}{\partial \timeVar} \phi(x, t) \d t \d x \\
			&\qquad+ \int\limits_{\supp{\phi}} \sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \frac{\partial}{\partial \spaceVar} \phi(x, t) \d t \d x \\
			&\qquad\qquad + \int\limits_{\supp{\phi}} \abs{\solution_{0}(\omega, x) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} \phi (x, 0) \d x .
		\end{align}
	\end{subequations}
	Taking the integral over $\supp{\phi}$ is a bounded linear operator, because $\supp{\phi}$ is compact.
	Since linear operators are bounded if and only if they are continuous, we know that taking integrals over a compact domain is a continuous operation.
	Now, since the composition of a continuous function with a measurable function is measurable by {\cite[Lemma 4.22]{Aliprantis2006Infinite}}, we obtain the measurability of all three integrals in $\omega \in \Omega$.
	Therefore, the measurability of $J^{\alpha}$ in $\omega \in \Omega$ follows and the proof is complete.
\end{proof}

We are now able to state the main result about the measurability of the solution map $u: \Omega \rightarrow L^{\infty}(\mathbb{R} \times \timeInterval) \cap C^{0}(\timeInterval; L^{1}_{\text{loc}}(\mathbb{R}))$.
\begin{theorem}[Measurability of the stochastic entropy solution]
	\label{thm:measurability}
	Let the flux function $\mathfrak{\solFlux}$ satisfy Assumptions \ref{as:stochFlux} and \ref{as:stochFlux_measurability} and let the initial condition satisfy $\solution_{0} \in L^p (\Omega, L^1(\mathbb{R}))$.
	Then, the mapping $u: \Omega \rightarrow L^{\infty}(\mathbb{R} \times \timeInterval) \cap C^{0}(\timeInterval; L^{1}_{\text{loc}}(\mathbb{R})), \omega \mapsto u(\omega, \cdot, \cdot)$ is strongly measurable, which means that it is measurable and $u (\Omega) \subset L^{\infty}(\mathbb{R} \times \timeInterval) \cap C^{0}(\timeInterval; L^{1}_{\text{loc}}(\mathbb{R}))$ is separable.
\end{theorem}

\begin{proof}
	For $N \in \mathbb{N}$, define the space
	\begin{equation}
		S_N \coloneqq \set{\phi \in C^{\infty}(\mathbb{R} \times \mathbb{T}; \mathbb{R}_{\geq0}): \supp(\phi) \subseteq \set{(x, t) \in \mathbb{R} \times \mathbb{T}: \abs{x}, t \leq N}} .
	\end{equation}
	Each space $S_N$ is a subspace of $C^{\infty}(\mathbb{R} \times \mathbb{T}; \mathbb{R}_{\geq0})$ and thus has a countable basis, since $C^{\infty}$ has a countable basis.
	Note that for every $\phi \in C^{\infty}_{c}(\mathbb{R} \times \mathbb{T}; \mathbb{R}_{\geq0})$, there exists $S_K$, $K \in \mathbb{N}$, such that $\phi \in S_K$.

	Now, for fixed $N \in \mathbb{N}$, let $(\varphi^{i}_{N}, i \in \mathbb{N}) \subset S_N$  be a basis of $S_N$.
	For fixed $i \in \mathbb{N}$ and $\alpha \in [M_0, \infty)$ (or $\alpha \in (-\infty, M_0]$), define the functional
	\begin{align*}
		J^{\alpha}_{i, N} &: \Omega \times \mathcal{X} \rightarrow \mathbb{R} , \\
		(\omega, \nu) &\mapsto \int \abs{\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} \frac{\partial}{\partial \timeVar} \varphi^{i}_{N} (x,t) \d t \d x \\
		&\hspace*{2cm} + \int \sign (\nu(x, t) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)) (\mathfrak{\solFlux}(\omega, x, \nu(x, t)) - \alpha) \frac{\partial}{\partial \spaceVar} \varphi^{i}_{N} (x,t) \d t \d x \\
		&\hspace*{3cm} + \int \abs{\solution_{0}(\omega, x) - \entropyConstant^{\pm}_{\steadyStatePar}(\stochVar, \spaceVar)} \varphi^{i}_{N} (x, 0) \d x .
	\end{align*}

	Due to Proposition \ref{prop:EntropyFunctional_Caratheodory}, the functional $J^{\alpha}_{i, N}$ is a Carathéodory map, i.e., it is measurable w.r.t. $\omega \in \Omega$ and continuous w.r.t. $w \in \mathcal{X}$.
	Note, this implies $-J^{\alpha}_{i, N}$ being a Carathéodory map.
	To show the measurability of $u$, we define the set-valued map $\Xi^{\alpha}_{i, N}: \Omega \rightrightarrows \mathcal{X}$
	\begin{align*}
		\Xi^{\alpha}_{i, N}(\omega) = \{ \nu \in \mathcal{X} | J^{\alpha}_{i, N}(\omega, \nu) \geq 0 \} , \text{ for } \omega \in \Omega .
	\end{align*}
	By {\cite[Prop. 6.3.4]{Geletu2006Introduction}}, this set-valued map $\Xi^{\alpha}_{i, N}(\omega)$ is measurable, since $\mathcal{X} = L^{\infty}(\mathbb{R} \times \timeInterval) \cap C^{0}(\timeInterval; L^{1}_{\text{loc}}(\mathbb{R}))$ is separable.

	Define now the set-valued map $\Xi^{\alpha}: \Omega \rightrightarrows \mathcal{X}$ via
	\begin{equation}
		\Xi^{\alpha}(\omega) = \bigcap_{i \in \mathbb{N}} \bigcap_{N \in \mathbb{N}} \Xi^{\alpha}_{i, N}(\omega) ,
	\end{equation}
	which is measurable as a countable intersection of measurable maps.

	We are now able to define the set-valued map $\Xi: \Omega \rightrightarrows \mathcal{X}$ as
	\begin{equation}
		\Xi(\omega) = \bigcap_{\alpha \in [M_0, \infty) \cap \mathbb{Q}} \Xi^{\alpha}(\omega) ,
	\end{equation}
	which is also measurable as a countable intersection of measurable maps.
	Note, if $\alpha \in (-\infty, M_0]$, we need to intersect over $\alpha \in (-\infty, M_0] \cap \mathbb{Q}$ instead of $\alpha \in [M_0, \infty) \cap \mathbb{Q}$.

	Note, since the adapted entropy functionals $J^{\alpha}_{i, N}$ depend continuously on $\alpha \in [M_0, \infty)$ it is sufficient to only consider $\alpha \in [M_0, \infty) \cap \mathbb{Q}$ here, since the rational numbers are dense in $\mathbb{R}$.
	To see this, let us first note that the intersection over $\alpha$ is necessary to select the adapted entropy solution, which satisfies the adapted entropy condition for all $\alpha \in [M_0, \infty)$.
	Let us now assume, that for $\alpha \in [M_0, \infty) \cap \mathbb{Q}$ the set $\Xi^{\alpha}(\omega)$ does contain the adapted entropy solution $u$ and an additional function $\tilde{u}$, satisfying the entropy condition for these $\alpha$.
	This implies, as a direct consequence of the existence and uniqueness theorem \ref{thm:stochExistenceUniqueness}, that for $\gamma \in [M_0, \infty) \setminus \mathbb{Q}$, the image of $\Xi^{\gamma}(\omega)$ does only contain the adapted entropy solution.
	By definition of the adapted entropy functional, this means that there exist $i_{-}, N_{-} \in \mathbb{N}$ such that
	\begin{align*}
		J^{\alpha}_{i, N}(\omega, \tilde{u}) &\geq 0 \qquad \forall i, N \in \mathbb{N} , \\
		J^{\gamma}_{i_{-}, N_{-}}(\omega, \tilde{u}) &< 0 \qquad \text{for } i_{-}, N_{-} \in \mathbb{N} .
	\end{align*}
	Since $J^{\gamma}_{i_{-}, N_{-}}$ is continuous in $\gamma$, we know that there exist an $\varepsilon$-neighbourhood $\mathcal{U}_{\varepsilon}$ of $J^{\gamma}_{i_{-}, N_{-}}(\omega, \tilde{u})$ such that for every $\xi \in \mathcal{U}_{\varepsilon}$, we have $\xi < 0$.
	From $\mathbb{Q}$ being dense in $\mathbb{R}$ and the continuity of $J^{\gamma}_{i_{-}, N_{-}}$ in $\gamma$, we know that there exist $\delta > 0$ and $\tilde{\gamma} \in [M_0, \infty) \cap \mathbb{Q}$, such that $\abs{\gamma - \tilde{\gamma}} < \delta$ and $J^{\tilde{\gamma}}_{i_{-}, N_{-}}(\omega, \tilde{u}) \in \mathcal{U}_{\varepsilon}$.
	However, $J^{\tilde{\gamma}}_{i_{-}, N_{-}}(\omega, \tilde{u}) \in \mathcal{U}_{\varepsilon}$ implies that $J^{\tilde{\gamma}}_{i_{-}, N_{-}}(\omega, \tilde{u}) < 0$.
	But this is a contradiction to $J^{\alpha}_{i, N}(\omega, \tilde{u}) \geq 0$ for all $i, N \in \mathbb{N}$ and all $\alpha \in [M_0, \infty) \cap \mathbb{Q}$.
	Consequently, it is sufficient to intersect over $\alpha \in [M_0, \infty) \cap \mathbb{Q}$ to select the adapted entropy condition.

	Due to the pathwise existence and uniqueness theorem \ref{thm:stochExistenceUniqueness}, we know that for fixed $\omega \in \Omega$, $\Xi(\omega)$ only contains the pathwise unique adapted entropy solution of \eqref{eq:stochasticConservationLaw}.
	Thus, the solution map $u: \Omega \rightarrow \mathcal{X}$ is measurable.
\end{proof}

\section{Stochastic jump coefficient}
\label{sec:stochJumpCoeff}

In this section, we focus on the specific choice of a multiplicative flux function (see also Example \ref{ex:multiplicativeFlux} and \ref{ex:multiplicativeFlux_measurability}), where the stochastic coefficient models heterogeneities in a medium.
For this jump coefficient, we use the random coefficient introduced in \cite{Barth2018study, Barth2020Numerical} for an elliptic diffusion problem.
It consists of a (spatial) Gaussian random field with additive discontinuities on random submanifolds.
Since the coefficient is intended to model fractured or heterogeneous media, it is assumed to be time-independent. \\
In general, it will not be possible to sample the coefficient exactly.
Therefore, we also introduce suitable approximations of the coefficient.

\begin{definition}[Stochastic jump coefficient]
	\label{def:stochasticJumpCoefficient}

	Let $\spaceDomain \subset \mathbb{R}$ be a compact subset.
	We define the stochastic jump coefficient as a function
	\begin{equation}
		\label{eq:stochasticJumpCoefficient}
		\begin{aligned}
			a: \Omega \times \mathbb{R} &\rightarrow \mathbb{R}_{> 0} \ , \\
			(\omega, x) &\mapsto \bar{\jumpCoeff}(x) + \phi (W(\omega, x)) + P (\omega, x) \ ,
		\end{aligned}
	\end{equation}
	where
	\begin{itemize}
		\item $\bar{\jumpCoeff} \in C(\mathbb{R}; \mathbb{R}_{\geq 0})$ is a deterministic, uniformly bounded mean function.
		\item $\phi \in C^{1}(\mathbb{R}; \mathbb{R}_{>0})$ is a continuously differentiable, positive mapping.
		\item For a (zero-mean) Gaussian random field $W_{\mathbb{R}} \in L^2(\Omega; L^2(\mathbb{R}))$ associated to a non-negative, symmetric trace class (covariance) operator $Q: L^2(\mathbb{R}) \rightarrow L^2(\mathbb{R})$, the truncated Gaussian random field $W \in L^{2}(\Omega; L^2(\mathbb{R}))$ is defined as
		\begin{equation}
			\label{eq:truncatedGaussField}
			\begin{aligned}
				W (\omega, x) =
				\begin{cases}
					W_{\mathbb{R}} (\omega, x) , & x \in \spaceDomain \\
					\min \left( W_{\mathbb{R}} (\omega, x), \sup_{x \in \spaceDomain} W_{\mathbb{R}} (\omega, x) \right) , & x \in \mathbb{R} \setminus \spaceDomain .
				\end{cases}
			\end{aligned}
		\end{equation}
		\item $\mathcal{T}: \Omega \rightarrow \mathcal{B}(\spaceDomain), \ \omega \mapsto \{ \mathcal{T}_1, \ldots, \mathcal{T}_{\tau} \}$ is a random partition of $\spaceDomain$, i.e., the $\mathcal{T}_i$ are disjoint open subsets of $\spaceDomain$ with $\overline{\spaceDomain} = \bigcup_{i=1}^{\tau} \overline{\mathcal{T}_i}$. The number of elements in $\mathcal{T}$ is a random variable $\tau: \Omega \rightarrow \mathbb{N}$ on $(\stochDomain, \sigmaAlgebra, \mathbb{P})$.\\
		For $\domainOfInterest_{l}$ and $\domainOfInterest_{r}$ being the left and right boundary of $\spaceDomain$, respectively, we define $\mathcal{T}_0 := (-\infty, \domainOfInterest_{l})$ and $\mathcal{T}_{\tau+1} := (\domainOfInterest_{r}, +\infty)$.
		\item A measure $\lambda$ on $\spaceDomain, \mathcal{B}(\spaceDomain)$ is associated to $\mathcal{T}$ and controls the position of the random elements $\mathcal{T}_i$.
		\item $(\jumpHeights{i}, i \in \mathbb{N}_{0})$ is a sequence of random variables on $(\stochDomain, \sigmaAlgebra, \mathbb{P})$ with arbitrary positive distribution(s) satisfying $\jumpHeights{i} < \infty$ for every $\omega \in \Omega$.
		Further, the sequence $(\jumpHeights{i}, i \in \mathbb{N}_{0})$ is independent of $\tau$ (but not necessarily i.i.d.) and we have
		\begin{equation}
			\label{eq:defJumpField}
			P: \Omega \times \spaceDomain \rightarrow \mathbb{R}_{>0}, \quad (\omega, x) \mapsto \sum_{i = 0}^{\tau+1} \boldsymbol{1}_{\mathcal{T}_i}(x) \jumpHeights{i}(\omega) \ .
		\end{equation}
	\end{itemize}
\end{definition}
\begin{remark}
	We note that we do not require the truncated Gaussian random field $W$ and the jump field $P$ to be stochastically independent.

	Further, the measure $\lambda$ associated to $\mathcal{T}$ does not only affect the average number of partition elements $\mathbb{E} (\tau)$, but also the size of the partition elements $\mathcal{T}_i$.
\end{remark}

\begin{corollary}[Assumptions on stochastic jump coefficient are satisfied]
	The stochastic jump coefficient defined in Definition \ref{def:stochasticJumpCoefficient} satisfies the assumptions for the stochastic flux function, see Example \ref{ex:multiplicativeFlux}.
\end{corollary}

\begin{proof}
	To satisfy the assumptions for a multiplicative flux function, the stochastic coefficient $a (\omega, x)$ has to be continuous in $x \in \mathbb{R}$ except on a closed set $\mathcal{N}$ of measure zero (which might depend on $\omega \in \Omega$).
	Further, the coefficient needs to have positive spatially bounded paths, where the bounds might be stochastic.

	Since the number of elements $\tau$ in the partition $\mathcal{T}$ is finite for every $\omega \in \Omega$, the jump field $P (\omega, \cdot)$ is a piecewise constant function with finitely many disconitnuities.
	Further, since $\spaceDomain$ is compact, the set $\mathfrak{D}(\omega)$ of discontinuities is closed and has measure zero.
	Thus, since the deterministic mean function $\bar{\jumpCoeff}$ is continuous by assumption and we can always consider a continuous modification of the truncated Gaussian random field $W$, it follows that $a (\omega, x)$ is continuous in $x \in \mathbb{R}$, except on a closed set $\mathfrak{D}(\omega, \cdot)$ of measure zero.

	It remains to show that the coefficient $a$ has positive spatially bounded paths:
	We have that $\jumpHeights{i} > 0$ for all $i \in \mathbb{N}_0$ by construction.
	Thus, there exists a $\jumpCoeff_{-} (\omega) > 0$, such that $P(\omega, x) \geq \jumpCoeff_{-} (\omega) > 0$ for every $x \in \mathbb{R}$.
	Since $\bar{\jumpCoeff}$ and $\phi$ are nonnegative, we have that $a(\omega, x) \geq \jumpCoeff_{-} (\omega) > 0$

	Further, since $\tau (\omega) < \infty$, for each $\omega \in \Omega$ there exists $P_{+} (\omega) < \infty$ with $P(\omega, x) \leq \max_{0 \leq i \leq \tau (\omega)+1} \jumpHeights{i}(\omega) \eqqcolon P_{+} (\omega)$.
	Since the Gaussian random field $W_{\mathbb{R}}(\omega, x)$ is continuous, it is bounded on $\spaceDomain$.
	Then, by construction, it follows that there exists $W_{+} < \infty$ such that $W (\omega, x) \leq W_{+}$.
	Since the mean function $\bar{\jumpCoeff}$ is bounded by assumption, we have shown that the stochastic jump coefficient from Definition \ref{def:stochasticJumpCoefficient} has positive spatially bounded paths, which proves the assertion.
\end{proof}
\begin{remark}[Alternative construction of the Gaussian random field]
	\label{rem:alternativeGaussFieldConstruction}
	The rather complicated construction of the random field $W$ via the domain of interest $\spaceDomain$ is necessary to ensure the boundedness of the spatial paths of the coefficient.
	Alternatively, it is possible to replace $\phi(W(\omega, x))$ by $\phi(W_{\mathbb{R}}(\omega, x))$ in Equation \eqref{eq:stochasticJumpCoefficient}, if one additionally assumes $\phi \in C^{1}(\mathbb{R}; \mathbb{R}_{>0})$ to be bounded.
\end{remark}

\subsection{Numerical approximation of the stochastic jump coefficient}

In general, the structure of the stochastic jump coefficient does not allow us to sample from its exact distribution.
To approximate the Gaussian field $W_{\mathbb{R}}$, one can use the Karhunen-Lo\`{e}ve expansion:\\
Let $((\eta_{i}, e_{i}), i \in \mathbb{N})$ denote a sequence of eigenpairs of the covariance operator $Q$.
Here, we assume the eigenvalues to be given in a decaying order, i.e., $\eta_{1} \geq \eta_{2} \geq \cdots \geq 0$.
As the covariance operator is trace class, the Gaussian random field $W_{\mathbb{R}}$ admits the representation
\begin{equation}
	\label{eq:KarhunenLoeve}
	W_{\mathbb{R}} = \sum_{i \in \mathbb{N}} \sqrt{\eta_{i}} e_{i} Z_{i} \ ,
\end{equation}
where $(Z, i \in \mathbb{N})$ is a sequence of independent standard normally distributed random variables.
An obvious approximation of the random field $W_{\mathbb{R}}$ is then given by the truncated Karhunen-Lo\`{e}ve expansion $\gaussField^{\KLindex}_{\mathbb{R}}$ of $W_{\mathbb{R}}$, which is defined as
\begin{equation}
	\label{eq:truncatedKarhunenLoeve}
	\gaussField^{\KLindex}_{\mathbb{R}} := \sum_{i = 1}^{N} \sqrt{\eta_{i}} e_{i} Z_{i} \ .
\end{equation}
Here, $N \in \mathbb{N}$ is called the \textit{cut-off index} of $\gaussField^{\KLindex}_{\mathbb{R}}$.
Now, with this approximation of $W_{\mathbb{R}}$, we can apply the truncation \eqref{eq:truncatedGaussField} to obtain an approximation of the truncated Gaussian random field $W$.
Note, if we use the construction proposed in Remark \ref{rem:alternativeGaussFieldConstruction} for $W$, instead of applying the truncation to $W_{\mathbb{R}}$, we can directly use the truncated Karhunen-Lo\'{e}ve expansion $\gaussField^{\KLindex}_{\mathbb{R}}$ as an approximation.

\section{Numerical Experiments}

In this section, we demonstrate, how the unknown $u$ may be approximated.
Therefore, in Section \ref{ssec:adaptiveDiscretization}, we introduce a sample-adapted Finite Volume scheme.
Here, sample-adaptivity means that the Finite Volume mesh is aligned \textit{a-priori} with the discontinuities of the stochastic jump coefficient $a$.
Consequently, the resulting discretization is stochastic, i.e., the mesh changes for each $\omega \in \Omega$.
This is in contrast to classical \textit{adaptive} Finite Volume methods, which are based on a-posteriori error estimates and remeshing the grid various times for each sample.

Afterwards, we investigate the approximations of the solution of the stochastic conservation law \eqref{eq:stochasticConservationLaw} with the specific flux function proposed in Example \ref{ex:multiplicativeFlux} combined with the stochastic jump coefficient introduced in Section \ref{sec:stochJumpCoeff}.
Therefore, we conduct experiments showing the performance of different numerical approximations and show the influence of different types of random fields on the pathwise solution of the stochastic conservation law.
We conclude this section by investigating the strong error of the solution and how it is affected by the various parameters of the stochastic jump coefficient.

To the best of the author's knowledge, writing a general numerical scheme to compute the adapted entropy solution in the sense of Audusse and Perthame \cite{Audusse2005Uniqueness} is still an open problem.
However, Godunov-type methods have been successfully applied to hyperbolic conservation laws with discontinuous flux functions for many other notions of solutions, see \cite{Adimurthi2005Godunov} for so-called solutions of type $(A, B)$ or \cite{Towers2000Convergence, Towers2018Convergence} for solutions satisfying an interface condition at the discontinuity points.
In 2019, Towers \cite{Towers2020existence} proved convergence of a finite difference scheme to the adapted entropy solution defined by Audusse \cite{Audusse2005Uniqueness} for a special class of conservation laws with flux functions in the space of bounded variation (BV).
Recently, Ghoshal et. al. \cite{Ghoshal2020Convergence} proved the convergence of a Godunov scheme to the Audusse-Perthame adapted entropy solution for the case of BV flux functions.

While the method proposed in \cite{Ghoshal2020Convergence} might be suitable to approximate the proposed stochastic conservation law for a broad class of stochastic jump coefficients, we note that in general, we cannot ensure our flux function to have bounded variation without restricting the possible choice of the covariance operator $Q$.

Throughout this section, we restrict ourselves to a bounded domain $\spaceDomain \subset \mathbb{R}$ and a finite time interval $\mathbb{T} = (0, T) \subset \mathbb{R}_{> 0}$.
Unless stated otherwise, we have $\spaceDomain \times \mathbb{T} = (0, 1)^2$.

\subsection{Sample-adapted discretization}
\label{ssec:adaptiveDiscretization}

In this section, we first propose different meshing strategies for approximating the adapted entropy solution via a Finite Volume scheme.
The main difference of these meshing strategies is, how they account for the discontinuities in the flux function.
Afterwards, we summarize the Finite Volume discretization for the reader's convenience.

\subsubsection{Sample-adapted meshing strategies}
\label{sssec:meshing-strategies}

In this section, we propose two meshing strategies that account for the discontinuities of the flux function.
First, we describe a samplewise jump-adapted meshing.
Such a meshing has already been used in \cite{Barth2020Numerical} for creating Finite Element meshes to approximate advection-diffusion problems.
Afterwards, we propose an improvement of this sample-adapted meshing for non-linear advection equations.
This novel strategy accounts for the standing-wave profiles occurring due to the flux discontinuities.

\paragraph{Samplewise jump-adapted meshing}
\label{par:jump-adaptive_meshing}

Let $\omega \in \Omega$ be arbitrary but fixed and denote by $\mathfrak{D} (\omega) \subset \mathbb{X}$ the set of discontinuities of $a(\omega, \cdot)$.
For $N_{\spaceVar} \in \mathbb{N}$, let $\spaceDomain_{\discretization} = \{ \gridPoint{i + \frac{1}{2}}\}_{i = 0}^{N_{\spaceVar}}$ be a decomposition of $\mathbb{X}$ with $\gridPoint{i - \frac{1}{2}} < \gridPoint{i+\frac{1}{2}}$ such that $\mathfrak{D} (\omega) \subset \spaceDomain_{\discretization}$.
We denote the $i$-th grid cell by $\gridCell{i} := (\gridPoint{i - \frac{1}{2}}, \gridPoint{i + \frac{1}{2}})$ and its associated spatial mesh size is defined as $\spaceStepSize{i} = \gridPoint{i+\frac{1}{2}} - \gridPoint{i - \frac{1}{2}}$.
Similarly, for $N_{\timeVar} \in \mathbb{N}$, let $\timeInterval_{\discretization} = \{ \timePoint{n} \}_{n = 0}^{N_{\timeVar}}$ be a discretization of the time interval with time step size $\timeStepSize{n} = \timePoint{n + 1} - \timePoint{n}$, for $n = 0, \ldots, N_{\timeVar} - 1$.
Note, the above construction ensures that each spatial discontinuity of the flux function is aligned with a cell interface of the Finite Volume method.

\paragraph{Samplewise jump-adapted wave-cell meshing}
\label{par:wave-cell_meshing}

The samplewise jump-adapted meshing method accounts for the discontinuities in the flux function, which produce standing wave profiles in the solution.
However, on meshes which allow for large step sizes, this standing wave profile might not be approximated sufficiently, if cells next to discontinuities are large.
Let again $\omega \in \Omega$ be arbitrary but fixed and let a samplewise jump-adapted mesh be given.
We denote by $\discontinuityInterface{i} \in \spaceDomain_{\discretization}$ the grid cell interfaces that correspond to a discontinuity of the flux function and by $\waveCell{i}^{l,r}$ its left and right grid cell.
Further, we denote by $\Delta \waveCell{i}$ the size of these discontinuity adjacent cells, which we will call \textit{wave cells}.
If
\begin{equation}
	\label{eq:waveCellCondition}
	\Delta \waveCell{i} \geq 2 \min_{i} \spaceStepSize{i} ,
\end{equation}
we propose to introduce another grid point, such that the new wave cell satisfies
\begin{equation}
	\label{eq:waveCellSizeCondition}
	\min_{i} \spaceStepSize{i} \leq \Delta \waveCell{i} < 2 \min_{i} \spaceStepSize{i} .
\end{equation}
As we see in the subsequent experiments, this meshing method can significantly improve the approximation of standing waves occurring in the solution.
Note, this new mesh does not change the Courant-Friedrichs-Lewy (CFL) condition of the jump-adapted mesh, as it does not affect the minimal step size of the discretization.

\subsubsection{Finite Volume method}
\label{sssec:FV-method}

In this section, we summarize the main ideas of the Finite Volume scheme as a spatial discretization.
The formulation is independent of the chosen meshing strategy of the previous section.
We employ a classical Finite Volume discretization in conservative form, which is given by
\begin{equation}
	\label{eq:FV_approximation}
	\approxSolution{n+1}{i} = \approxSolution{n}{i} - \frac{\timeStepSize{n}}{\spaceStepSize{i}} \left( \outFlux{}{i} - \inFlux{}{i} \right) \ , \qquad
	\approxSolution{0}{i} \approx \frac{1}{\spaceStepSize{i}} \int_{\gridCell{i}} \solution_{0}(\omega, x) \d{x} \ .
\end{equation}

In the presented numerical experiments, we employ the Godunov flux for spatially dependent flux functions, which is given by
\begin{equation}
	\label{eq:GeneralizedGodunovFlux}
	\inFlux{}{i} =
	\begin{cases}
		\min\limits_{\approxSolution{}{i-1} \leq \theta \leq \approxSolution{}{i}} \mathfrak{\solFlux} (\omega, \gridPoint{i - \frac{1}{2}}, \theta) \qquad \text{if } \approxSolution{}{i-1} \leq \approxSolution{}{i} \ , \\
		\max\limits_{\approxSolution{}{i-1} \leq \theta \leq \approxSolution{}{i}} \mathfrak{\solFlux} (\omega, \gridPoint{i - \frac{1}{2}}, \theta) \qquad \text{if } \approxSolution{}{i} \leq \approxSolution{}{i-1} \ .
	\end{cases}
\end{equation}

Note, in the case of a multiplicative flux function $ \mathfrak{\solFlux} (\omega, x, u) = a (\omega, x) f(u)$, with $f$ being convex, we can simplify Equation \eqref{eq:GeneralizedGodunovFlux}:
\begin{equation}
	\label{eq:unimodalGodunovFlux_convex}
	\inFlux{}{i} = \max \left\{ a(\omega, \gridPoint{i-1}) f(\max(\approxSolution{}{i-1}, 0)), a(\omega, \gridPoint{i}) f(\min(\approxSolution{}{i}, 0)) \right\}.
\end{equation}
Note, when $a(\omega, \gridPoint{i-1}) = a (\omega, \gridPoint{i})$, Equation \eqref{eq:unimodalGodunovFlux_convex} reduces to the classical Godunov flux, which is used for conservation laws, where the flux function is not spatially dependent.

For the numerical time integration, we use the forward Euler time integration scheme, resulting the time marching equation
\begin{equation}
	\label{eq:timeMarching}
	\approxSolution{n+1}{i} = \approxSolution{n}{i} - \frac{\timeStepSize{n}}{\spaceStepSize{i}} \left( \outFlux{n}{i} - \inFlux{n}{i} \right)
\end{equation}
For stability reasons, we assume that the discretizations $\spaceDomain_{\discretization}$ and $\timeInterval_{\discretization}$ satisfy the Courant-Friedrichs-Lewy stability condition:
\begin{equation}
	\label{eq:CFLcondition}
	\frac{\timeStepSize{n}}{\min\limits_{i} \spaceStepSize{i}} L_u^n < 1 , \qquad \text{for all } n = 0, \ldots, N_{\timeVar} - 1
\end{equation}
where
\begin{equation}
	\label{eq:CFLconstant_general}
	L_u^n = \max_{i} \set{\abs{\frac{\partial \mathfrak{\solFlux}(\omega, \gridPoint{i - \frac{1}{2}}, \approxSolution{n}{i})}{\partial \approxSolution{}{}}}}
\end{equation}
in the general case, and
\begin{equation}
	\label{eq:CFLconstant_multiplicative}
	L_u^n = \norm{a (\omega)}_{\infty} \cdot \max_{i} \set{\abs{f'(\approxSolution{n}{i})}}
\end{equation}
in case of a multiplicative flux function.

\subsection{Stochastic Burgers' equation}

For the numerical experiments described in this section, we consider the stochastic Burgers' equation, which is given by
\begin{equation}
	\label{eq:stochasticBurgers}
	\begin{aligned}
		u_{t} + \left( a(\omega, x) \frac{u^2}{2} \right)_{x} &= 0 &&\text{in } \Omega \times \mathbb{X} \times \mathbb{T} \ , \\
		u(\omega, x, 0) &= \solution_{0} (\omega, x) &&\text{in } \Omega \times \mathbb{X} \times \set{0} \ ,
	\end{aligned}
\end{equation}
which is equipped with periodic boundary conditions, i.e., $a (\omega, 0) \frac{u (\omega, 0, t)^2}{2} = a (\omega, 1) \frac{u (\omega, 1, t)^2}{2}$, unless stated otherwise.

The deterministic mean value of the stochastic jump coefficient is given by $\bar{\jumpCoeff} \equiv 0$ and we choose $\phi (\xi) = \exp(\xi)$.
The Gaussian field $W_{\mathbb{R}}$ is characterized by the \textit{Mat\'{e}rn covariance operator} with smoothness parameter $\nu > 0$, variance $\sigma^2 > 0$ and correlation length $\rho > 0$, given by
\begin{equation}
	\label{eq:MaternCovariance}
	\begin{aligned}
		Q_{M}&: L^2(\spaceDomain) \rightarrow L^2(\spaceDomain) \ ,\\
		[Q_{M} \varphi] (y) &:= \int_{\mathbb{X}} \sigma^2 \frac{2^{1-\nu}}{\Gamma(\nu)} \left( \sqrt{2 \nu} \frac{\abs{x - y}}{\rho} \right)^{\nu} K_{\smoothness} \left( \sqrt{2 \nu} \frac{\abs{x - y}}{\rho} \right) \varphi(x) dx \quad \forall \varphi \in L^2(\spaceDomain) \ .
	\end{aligned}
\end{equation}
Here, $\Gamma$ denotes the Gamma function and $K_{\smoothness}$ denotes the modified Bessel function of second kind with $\nu$ degress of freedom.
The spectral basis of the covariance operator $Q_{M}$ may be approximated via Nystr\"{o}m's method \cite{Williams2006Gaussian}.

\subsection{Parameter study}

In this section, we investigate the influence of various parameters of the stochastic jump coefficient on the solution of Equation \eqref{eq:stochasticBurgers} and on the convergence rate of its approximation.
Therefore, we consider the strong $L^{1}$-error, i.e.,
\begin{equation}
	\label{eq:strongError_L1}
	\mathcal{E} (\omega) \coloneqq \norm{u^{\text{ref}}_{\Delta} (\omega, \cdot, T) - u_{\Delta} (\omega, \cdot, T)}_{L^{1}(\mathbb{X})} .
\end{equation}
Here, $u^{\text{ref}}_{\Delta}$ is a numerical reference solution computed on a finer grid and $u_{\Delta}$ is the considered approximation.
We first investigate the influence of the jump field.
Afterwards, we present two experiments providing insight about the performance of the explicit and implicit time integration scheme for various simulation settings.
We conclude by showing the influence of the parameters of the Gaussian random field on the numerical approximations. \\
Note, even though we only consider the $L^{1}$-error in this parameter study, the findings qualitatively also remain true, when we consider the $L^{2}$-error instead.

\subsubsection{Jump field parameters}
\label{sssec:jumpFieldParameters}

\paragraph{Distance between jumps}
In this experiment, we investigate, how the distance between two jumps influences the performance of the finite volume approximation.
Here, we also compare the sample-adapted discretization with the corresponding equidistant approximation.
Therefore, we denote the jump area with size $\delta > 0$ as $J_{\delta} = (1 - \frac{\pi}{10} - \frac{\delta}{2}, 1 - \frac{\pi}{10} + \frac{\delta}{2}) \subset (0, 1)$ and define our jump coefficient as $a(x) = P(x)$, where $P(x)$ has one of the following forms:
\begin{equation}
	P^{\text{up}}(x) = \begin{cases} \frac{1}{2} & \text{for } x \in (0, 1) \setminus J_{\delta} \\ \frac{3}{2\delta} & \text{for } x \in J_{\delta} \end{cases} , \qquad P^{\text{down}}(x) = \begin{cases} \frac{3}{2} & \text{for } x \in (0, 1) \setminus J_{\delta} \\ \delta & \text{for } x \in J_{\delta} \end{cases}.
\end{equation}
In Figures \ref{fig:jumpDistance_jumpUp} and \ref{fig:jumpDistance_jumpDown}, the convergence rates of the finite volume approximation is plotted for $P^{\text{up}}$ and $P^{\text{down}}$, respectively, and for various sizes of the jump area $J_{\delta}$.
\begin{figure}[htbp!]
	\begin{subfigure}{.32\textwidth}
		\includegraphics{./pictures/jumpInfluence/jumpDistance_jump-up_L1/delta_2-4.pdf}
		\caption{$\delta = 2^{-4}$.}
	\end{subfigure}
	\hfill
	\begin{subfigure}{.32\textwidth}
		\includegraphics{./pictures/jumpInfluence/jumpDistance_jump-up_L1/delta_2-6.pdf}
		\caption{$\delta = 2^{-6}$.}
	\end{subfigure}
	\hfill
	\begin{subfigure}{.32\textwidth}
		\includegraphics{./pictures/jumpInfluence/jumpDistance_jump-up_L1/delta_2-8.pdf}
		\caption{$\delta = 2^{-8}$.}
	\end{subfigure}
	\caption{Convergence rates for $a(x) = P^{\text{up}}(x)$ and various $\delta$.}
	\label{fig:jumpDistance_jumpUp}
\end{figure}
\begin{figure}[htbp!]
	\begin{subfigure}{.32\textwidth}
		\includegraphics{./pictures/jumpInfluence/jumpDistance_jump-down_L1/delta_2-4.pdf}
		\caption{$\delta = 2^{-4}$.}
	\end{subfigure}
	\hfill
	\begin{subfigure}{.32\textwidth}
		\includegraphics{./pictures/jumpInfluence/jumpDistance_jump-down_L1/delta_2-6.pdf}
		\caption{$\delta = 2^{-6}$.}
	\end{subfigure}
	\hfill
	\begin{subfigure}{.32\textwidth}
		\includegraphics{./pictures/jumpInfluence/jumpDistance_jump-down_L1/delta_2-8.pdf}
		\caption{$\delta = 2^{-8}$.}
	\end{subfigure}
	\caption{Convergence rates for $a(x) = P^{\text{down}}(x)$ and various $\delta$.}
	\label{fig:jumpDistance_jumpDown}
\end{figure}

While the convergence rate of the jump-adapted finite volume discretization does not seem to be affected by the distance between two jumps, the equidistant finite volume method has problems with approximating small distances.

\paragraph{Number of jumps}
We continue our numerical experiments by investigating, how the number of jumps in the coefficient affects the convergent rate of the (jump-adapted) finite volume discretization.
Therefore, we consider again the case of $a(x) = P(x)$, where $P(x)$ is defined as in Equation \eqref{eq:defJumpField},
\begin{equation}
	P: \spaceDomain \rightarrow \mathbb{R}_{>0}, \quad x \mapsto \sum_{i = 0}^{\tau+1} \boldsymbol{1}_{\mathcal{T}_i}(x) \jumpHeights{i} \ , \qquad P_i = \begin{cases} \frac{1}{2} & i \text{ odd} \\ \frac{3}{2} & i \text{ even} \end{cases} .
\end{equation}
Above, $\tau \in \mathbb{N}$ is the number of jumps and we set the jump positions $\jumpPosition{i}$ to be uniformly distributed over the domain, i.e., $\jumpPosition{i} \sim \mathcal{U}((0, 1))$.
Figure \ref{fig:numberOfJumps} shows the convergence rates for various values of $\tau$.
\begin{figure}[htbp!]
	\begin{subfigure}{.32\textwidth}
		\includegraphics{./pictures/jumpInfluence/numberOfJumps_L1/numJumps-4}
		\caption{$\tau = 4$ jumps.}
	\end{subfigure}
	\hfill
	\begin{subfigure}{.32\textwidth}
				\includegraphics{./pictures/jumpInfluence/numberOfJumps_L1/numJumps-16}
		\caption{$\tau = 16$ jumps.}
	\end{subfigure}
	\hfill
	\begin{subfigure}{.32\textwidth}
				\includegraphics{./pictures/jumpInfluence/numberOfJumps_L1/numJumps-64}
		\caption{$\tau = 64$ jumps.}
		\label{sfig:64jumps}
	\end{subfigure}
	\caption{Convergence rates for $a(x) = P(x)$ with varying number of jumps.}
	\label{fig:numberOfJumps}
\end{figure}

One can see that an increasing number of jumps leads to a slower convergence rate.
Another result is that the jump-adapted finite volume approximation does not improve the convergence rate of the method.
However, for a coefficient with many jumps, Figure \ref{sfig:64jumps} shows that the jump-adapted discretization might lead to a better constant.

\subsubsection{Explicit vs. implicit time integration}

In this section, we conduct experiments on the performance of the explicit and implicit Euler scheme to integrate the jump-adapted finite volume discretization.
Specifically, we investigate the influence of two parameters on the performance of the methods, measured in time-to-error:
(i) The influence of $\frac{\max (\Delta x)}{\min (\Delta x)}$, which influences the possible time step size due to the CFL condition in the explicit scheme and the lack thereof in the implicit scheme.
(ii) The influence of the estimated wave speed of the solution, which also influences the possible time step size.
In both experiments, we compare the behaviour for multiple discretization parameters.

\paragraph{Influence of $\frac{\max (\Delta x)}{\min (\Delta x)}$}

In this experiment, we consider the following coefficient for some parameter $\delta > 0$
\begin{equation}
	a(x) = \begin{cases}
		0.5 & \text{for } 0 < x < \frac{\pi}{5} - \frac{\delta}{2} \\
		1.5 & \text{for } \frac{\pi}{5} - \frac{\delta}{2} < x < \frac{\pi}{5} + \frac{\delta}{2} \\
		0.5 & \text{for } \frac{\pi}{5} + \frac{\delta}{2} < x < 1
	\end{cases} .
\end{equation}
Above, $\delta$ denotes the distance between the two jumps of $a$ and thus imposes a restriction on $\min (\Delta x)$.
The corresponding time-to-error plots for various values of $\delta$ are shown in Figure \ref{fig:explciitVsImplicit_DeltaX}.
\begin{figure}[htbp!]
	\begin{subfigure}{.325\textwidth}
		\includegraphics{./pictures/explicitVsImplicit/minStepSize_L1/delta_1e-2}
		\caption{$\delta = 10^{-2}$.}
	\end{subfigure}
	\hfill
	\begin{subfigure}{.325\textwidth}
		\includegraphics{./pictures/explicitVsImplicit/minStepSize_L1/delta_1e-3}
		\caption{$\delta = 10^{-3}$.}
	\end{subfigure}
	\hfill
	\begin{subfigure}{.325\textwidth}
		\includegraphics{./pictures/explicitVsImplicit/minStepSize_L1/delta_1e-4}
		\caption{$\delta = 10^{-4}$.}
	\end{subfigure}
	\caption{Time to error plot of the forward Euler and backward Euler time discretization of the jump-adapted finite volume approximation for different jump distances $\delta > 0$.}
	\label{fig:explciitVsImplicit_DeltaX}
\end{figure}

It is not surprising that the backward Euler time discretization needs more time with decreasing error, since the size of the nonlinear system, which needs to be solved at each time step, increases.
Note that the computation time does not increase with finer discretizations.
However, the amount of time necessary to evolve the solution over time with the explicit scheme depends highly on the allowed maximum time step size, which is directly influenced by $\delta$ through the CFL condition.

\paragraph{Influence of the estimated wave speed}

In this experiment, we conduct experiments on the influence of the estimated wave speed on the performance of the explicit and implicit time integrator, respectively.
Therefore, we consider an initial condition $\solution_{0} (x) = \kappa \sin (\pi x)$ with $\kappa \in \mathbb{R}_{>0}$ and a jump coefficient of the form
\begin{equation}
	a(x) = \begin{cases}
		P_1 & \text{for } 0 < x < \frac{\pi}{5} - \frac{1}{20} \\
		P_2 & \text{for } \frac{\pi}{5} - \frac{1}{20} < x < \frac{\pi}{5} + \frac{1}{20} \\
		P_1 & \text{for } \frac{\pi}{5} + \frac{1}{20} < x < 1
	\end{cases} .
\end{equation}
Figure \ref{fig:explicitVsImplicit_waveSpeed} visualizes the time-to-error plots for the explicit and implicit time integrator for different values of $\kappa$ and $P_1, P_2$.
\begin{figure}[h!tbp]
	\centering
	\begin{subfigure}[b]{.325\textwidth}
		\includegraphics{./pictures/explicitVsImplicit/waveSpeed_L1/timeToError_maxCoeff-20}
		\caption{$\kappa = 0.3$, $P_1 = 10.5$, $P_2 = 20$.}
	\end{subfigure}
	\hfill
	\begin{subfigure}[b]{.325\textwidth}
		\includegraphics{./pictures/explicitVsImplicit/waveSpeed_L1/timeToError_maxCoeff-50}
		\caption{$\kappa = 0.3$, $P_1 = 1$, $P_2 = 50$.}
	\end{subfigure}
	\hfill
	\begin{subfigure}[b]{.325\textwidth}
		\includegraphics{./pictures/explicitVsImplicit/waveSpeed_L1/timeToError_maxCoeff-50_sin-1}
		\caption{$\kappa = 1$, $P_1 = 1$, $P_2 = 50$.}
		\label{sfig:waveSpeed_50_sin1}
	\end{subfigure}
	\caption{Time to error plots of the forward and backward Euler scheme with a jump-adapted finite volume discretization for different values of $\kappa$ and $P_1, P_2$.}
	\label{fig:explicitVsImplicit_waveSpeed}
\end{figure}

As one can see in Figure \ref{sfig:waveSpeed_50_sin1}, the implicit time integration only outperforms the explicit time integration in the case of a coarse spatial discretization combined with a high maximum value of the jump coefficient.

\subsubsection{Parameters of Gaussian random field}

In this section, we investigate the influence of the parameters of the Gaussian random field on the convergence rates of the equidistant finite volume approximation.
Therefore, we set $a(\omega, x) = \exp (W_{\mathbb{R}}(\omega, x))$, i.e., we consider a log-Gaussian random field.
The Gaussian random field $W_{\mathbb{R}}$ is characterized by the Mat{\'e}rn covariance operator (see Equation \eqref{eq:MaternCovariance}) with smoothness parameter $\nu > 0$, variance $\sigma^2 > 0$ and correlation length $\rho > 0$.
In the subsequent experiments, we focus on the influence of the smoothness $\nu$ and the correlation length $\rho$. \\
Since the Gaussian field is stochastic, the Figures \ref{fig:Gaussian_smoothness} and \ref{fig:Gaussian_correlationLength} show convergence rates that are averaged over $20$ samples.
For both experiments, the variance was fixed at $\sigma^2 = 0.1$.
While the smoothness parameter $\nu$ was varied in the first experiment, the correlation length was fixed at $\rho = 0.1$ and in the second experiment, we set $\nu = \infty$ and varied the correlation length.
\begin{figure}[h!tbp]
	\centering
	\begin{subfigure}[b]{.32\textwidth}
		\includegraphics{./pictures/GaussianInfluence/smoothness_L1/smoothness-0-5}
		\caption{Mat{\'e}rn smoothness $\nu = 0.5$.}
	\end{subfigure}
	\hfill
	\begin{subfigure}[b]{.32\textwidth}
		\includegraphics{./pictures/GaussianInfluence/smoothness_L1/smoothness-1}
		\caption{Mat{\'e}rn smoothness $\nu = 1$.}
	\end{subfigure}
	\hfill
	\begin{subfigure}[b]{.32\textwidth}
		\includegraphics{./pictures/GaussianInfluence/smoothness_L1/smoothness-infty}
		\caption{Mat{\'e}rn smoothness $\nu = \infty$.}
	\end{subfigure}
	\caption{Convergence of the finite volume approximation for different smoothness parameters $\nu$ of the Mat{\'e}rn covariance kernel. The error estimation is based on $20$ samples.}
	\label{fig:Gaussian_smoothness}
\end{figure}

As one can see, with decreasing smoothness parameter, the convergence rate of the finite volume approximation gets lower.
This is, what one would expect, since rough paths of the Gaussian random field are not resolved as good as smooth paths.
A similar effect can be seen by reducing the correlation length of the covariance kernel, which is shown in Figure \ref{fig:Gaussian_correlationLength}.
\begin{figure}[h!tbp]
	\centering
	\begin{subfigure}[b]{.32\textwidth}
		\includegraphics{./pictures/GaussianInfluence/correlationLength_L1/correlationLength-0-01_L1}
		\caption{Correlation length $\rho = 0.01$.}
	\end{subfigure}
	\hfill
	\begin{subfigure}[b]{.32\textwidth}
		\includegraphics{./pictures/GaussianInfluence/correlationLength_L1/correlationLength-0-05_L1}
		\caption{Correlation length $\rho = 0.05$.}
	\end{subfigure}
	\hfill
	\begin{subfigure}[b]{.32\textwidth}
		\includegraphics{./pictures/GaussianInfluence/correlationLength_L1/correlationLength-0-1_L1}
		\caption{Correlation length $\rho = 0.1$.}
	\end{subfigure}
	\caption{Convergence of the finite volume approximation for different correlation length $\rho$ of the Mat{\'e}rn covariance kernel. The error estimation is based on $20$ samples.}
	\label{fig:Gaussian_correlationLength}
\end{figure}

\subsection{Pathwise convergence study}

In this section, we investigate the strong error of the jump-adapted and wave cell finite volume method when computing the stochastic Burgers equation \eqref{eq:stochasticBurgers}.
In the experiments described in Sections \ref{sssec:alternatingJumpField} and \ref{sssec:PoissonJumpField}, we consider specific jump coefficients, which have the form described in \eqref{eq:stochasticJumpCoefficient}.
Afterwards, in Section \ref{sssec:InclusionField}, we present experiments demonstrating the ability of the wave cell finite volume method.
Therefore, we consider a stochastic jump coefficient without the Gaussian random field part, i.e., having the simpler form $a (\omega, x) = \bar{\jumpCoeff} (x) + P (\omega, x)$.
For the convergence study, we consider both the strong $L^{1}$-error and the strong $L^{2}$-error, i.e.,
\begin{equation}
	\label{eq:strongError}
	\mathcal{E} (\omega) \coloneqq \norm{u^{\text{ref}}_{\Delta} (\omega, \cdot, T) - u_{\Delta} (\omega, \cdot, T)}_{L^{\star}(\mathbb{X})} .
\end{equation}
Here, $u^{\text{ref}}_{\Delta}$ is a numerical reference solution computed on a finer grid, $u_{\Delta}$ is the considered approximation and $\norm{\cdot}_{L^{\star}(\mathbb{X})}$ denotes either the $L^{1}(\mathbb{X})$ or the $L^{2}(\mathbb{X})$ norm.

\subsubsection{Alternating jump field with exponential Gaussian random field}
\label{sssec:alternatingJumpField}

In this first experiment, we consider a Gaussian random field $W$ with the Matérn covariance operator \eqref{eq:MaternCovariance} with a smoothness parameter $\nu = \frac{1}{2}$.
The jump field $P (\omega, x)$ has $\tau \sim Poi(5)+1$ jumps, which are located at the jump locations $\jumpPosition{i} \sim \mathcal{U} (\mathcal{D})$.
The jump heights $\jumpHeights{i}$ corresponding to the partition, are given by
\begin{align}
	\label{eq:alternatingJumpHeights}
	\jumpHeights{i} \sim
	\begin{cases}
		\mathcal{U} ([ \frac{1}{4}, \frac{3}{4}]) &\text{for } i \text{ odd}, \\
		\mathcal{U} ([ \frac{5}{4}, \frac{7}{4}]) &\text{for } i \text{ even}.
	\end{cases}
\end{align}
To illustrate this coefficient, Figure \ref{fig:alternatingExponential_samples} shows two realizations of the coefficient.
One can see that the occurring jumps are relatively small compared to the variation of the Gaussian random field, which is influenced by the small smoothness parameter $\nu$, the variance $\sigma^2$ and the correlation length $\rho$.
\begin{figure}[h!tbp]
	\centering
	\begin{subfigure}[b]{0.495\textwidth}
		\centering
		\includegraphics{./pictures/coefficientSamples/alternatingExponential/alternatingExponential-Sample-1}
	\end{subfigure}
	\hfill
	\begin{subfigure}[b]{0.495\textwidth}
		\centering
		\includegraphics{./pictures/coefficientSamples/alternatingExponential/alternatingExponential-Sample-2}
	\end{subfigure}
	\caption{Two samples of the alternating jump field with an exponential Gaussian random field.}
	\label{fig:alternatingExponential_samples}
\end{figure}

In Figure \ref{fig:alternatingExponential_errorVsStepSize}, the pathwise $L^{1}$- and $L^{2}$-error are shown for $50$ samples.
The corresponding reference solutions were computed with the wave cell finite volume method using a grid with approximately $4$ times as many discretization points.
One can see that the difference of the methods in not significant and that the order of convergence is similar.
Since the distance between the jumps is quiet large on average, the discontinuities are already resolved on relatively coarse discretizations.
Thus, the different in the discretizations has no significant effect and the presented behaviour is exactly, what we would expect.
\begin{figure}[h!tbp]
	\centering
	\begin{subfigure}[b]{.49\textwidth}
		\centering
		\includegraphics{./pictures/convergenceResults/alternatingExponential/L1-error_vs_invStepSize}
		\caption{Pathwise $L^{1}$ error.}
	\end{subfigure}
	\hfill
	\begin{subfigure}[b]{.49\textwidth}
		\centering
		\includegraphics{./pictures/convergenceResults/alternatingExponential/L2-error_vs_invStepSize}
		\caption{Pathwise $L^{2}$ error.}
	\end{subfigure}
	\caption{Strong error for the alternating jump field with exponential Gaussian random field. The error estimation is based on $50$ samples.}
	\label{fig:alternatingExponential_errorVsStepSize}
\end{figure}

\subsubsection{Poisson jump field with squared-exponential Gaussian random field}
\label{sssec:PoissonJumpField}

In this experiment, we consider a similar jump field as in Section \ref{sssec:alternatingJumpField}.
The considered jump field $P$ has $\tau \sim Poi(5)+1$ jumps at jump positions $\jumpPosition{i} \sim \mathcal{U} (\mathcal{D})$.
The corresponding jump heights are $\jumpHeights{i} \sim Poi(5) + 1$ distributed.
This leads to a significant difference in the jump field, since the jump heights are not globally bounded.
We additively combine this random jump field with a squared-exponential Gaussian random field, i.e., the Matérn covariance kernel \eqref{eq:MaternCovariance} having a smoothness parameter of $\nu = \infty$.

Two samples of this coefficient are shown in Figure \ref{fig:PoissonSquaredExponential_samples}.
Note that, in contrast to the previous random field, the jumps are higher and have a more significant effect on the sample compared to the variation of the Gaussian part.
Also, let us mention that this coefficient is more likely to have larger values than the random field presented in the previous subsection.
\begin{figure}[h!tbp]
	\centering
	\begin{subfigure}[b]{0.495\textwidth}
		\centering
		\includegraphics{./pictures/coefficientSamples/PoissonSquaredExponential/PoissonSquaredExponential-Sample-1}
	\end{subfigure}
	\hfill
	\begin{subfigure}[b]{0.495\textwidth}
		\centering
		\includegraphics{./pictures/coefficientSamples/PoissonSquaredExponential/PoissonSquaredExponential-Sample-2}
	\end{subfigure}
	\caption{Two realizations of the Poisson jump field with squared-exponential Gaussian random field.}
	\label{fig:PoissonSquaredExponential_samples}
\end{figure}

The pathwise $L^{1}$- and $L^{2}$-error are shown in Figure \ref{fig:PoissonSquaredExponential}.
The convergence behaviour and the absolute error are again very similar, which is as expected, since the jumps are resolved on relatively coarse grids.
However, we note a slightly higher convergence rate for both adapted discretization methods.
This is not surprising, since the heights of the standing shock profiles resulting from the discontinuities in the flux function are better approximated with the adapted methods.
\begin{figure}[h!tbp]
	\centering
	\begin{subfigure}[b]{.49\textwidth}
		\centering
		\includegraphics{./pictures/convergenceResults/PoissonSquaredExponential/L1-error_vs_invStepSize}
		\caption{Pathwise $L^{1}$ error.}
	\end{subfigure}
	\hfill
	\begin{subfigure}[b]{.49\textwidth}
		\centering
		\includegraphics{./pictures/convergenceResults/PoissonSquaredExponential/L2-error_vs_invStepSize}
		\caption{Pathwise $L^{2}$ error.}

	\end{subfigure}
	\caption{Strong error for the Poisson jump field with squared-exponential Gaussian random field estimated with $50$ samples.}
	\label{fig:PoissonSquaredExponential}
\end{figure}

\subsubsection{Jump field with random inclusions}
\label{sssec:InclusionField}

In this last experiment, we demonstrate the ability of the wave cell finite volume discretization.
Therefore, we consider the following stochastic jump coefficient having random inclusions:
Let $\tau \sim Poi(10)+1$ denote the number of inclusions with positions $\chi_i \sim \mathcal{U} (\mathcal{D})$.
The corresponding inclusion sizes are given by $\delta_i \sim \mathcal{U} ([10^{-5}, 10^{-3}])$ and the heights of the inclusions are given by
\begin{align}
	\label{eq:inclusionHeights}
	H_k \sim
	\begin{cases}
		Poi(30) &\text{for } k = 0, \\
		\frac{1}{\xi}, \xi \sim Poi(30) &\text{for } k = 1.
	\end{cases}
\end{align}
Let now $\mathfrak{I} \subset \mathcal{D}$ denote the union of all inclusion.
We define the jump field with random inclusions as
\begin{align}
	\label{eq:inclusionField}
	a (\omega, x) =
	\begin{cases}
		H_k &\text{for } x \in \mathfrak{I}, \\
		1 &\text{otherwise,}
	\end{cases}
\end{align}
where $k \sim B(1, \frac{1}{2})$ is a Bernoulli-$\frac{1}{2}$-distributed random variable, defining whether the inclusion height is given by $H_1$ or by $H_2$.

The random field is illustrated in Figure \ref{fig:inclusionField_samples}, where two samples of the coefficient are shown.
Note, the inclusions are marked with a circle to increase the visibility.
Also, we see that the character of the random field is very different to the previous ones, since no Gaussian field is involved and thus the coefficient does not vary away from the jumps.
However, due to the small size of the inclusions, this coefficient is promising to demonstrate the qualitative difference between the proposed meshing methods.
\begin{figure}[h!tbp]
	\centering
	\begin{subfigure}[b]{0.495\textwidth}
		\centering
		\includegraphics{./pictures/coefficientSamples/inclusionField/inclusionField-Sample-1}
	\end{subfigure}
	\hfill
	\begin{subfigure}[b]{0.495\textwidth}
		\centering
		\includegraphics{./pictures/coefficientSamples/inclusionField/inclusionField-Sample-2}
	\end{subfigure}
	\caption{Two realizations of the jump field with random inclusions.}
	\label{fig:inclusionField_samples}
\end{figure}

In Figure \ref{fig:inclusionField}, we can see the performance of the different discretization methods measured by the strong $L^{1}$ and $L^{2}$ error for $100$ samples.
While the adapted discretizations are converging with a similar rate, the standard finite volume method does not seem to converge at all or at a very low convergence rate.
This behaviour is what one would expect for this jump coefficient and the standard finite volume method:
Due to the small size of the inclusions, the standard finite volume scheme does not resolve the discontinuities in the flux function.
Thus, when refining the spatial step size, the standard finite volume discretization needs a much finer resolution of the grid to resolve the discontinuities than the adaptive discretization schemes.
Note, asymptotically, we expect all three discretization schemes to converge with the same convergence rate.  \\
We also note the better error constant in the wave cell finite volume method compared to the jump-adapted finite volume method.
This behaviour can be explained by the better approximation of the standing wave profiles resulting from the discontinuities of the flux function.
\begin{figure}[h!tbp]
	\centering
	\begin{subfigure}[b]{.49\textwidth}
		\includegraphics{./pictures/convergenceResults/inclusionField/L1-error_vs_invStepSize}
		\caption{Pathwise $L^{1}$ error.}
	\end{subfigure}
	\hfill
	\begin{subfigure}[b]{.49\textwidth}
		\includegraphics{./pictures/convergenceResults/inclusionField/L2-error_vs_invStepSize}
		\caption{Pathwise $L^{2}$ error.}
	\end{subfigure}
	\caption{Strong error for the stochastic jump field with random inclusions based on $100$ samples.}
	\label{fig:inclusionField}
\end{figure}

\section*{Declaration of interest}
None.

\section*{Acknowledgement}
The research leading to these results has received funding from the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) under Germany's Excellence Strategy EXC-2075 390740016 at the University of Stuttgart and it is gratefully acknowledged.

\printbibliography
\end{document}