1 Introduction
The usual hydrological model for contamination/tracer transport through a porous media is given by a second order advection-dispersion equation (ADE) of the form $\frac{\partial c}{\partial t}=-\frac{\partial }{\partial x}[v(x)c]+\frac{{\partial ^{2}}}{\partial {x^{2}}}[D(x)c]$, where $c(x,t)$ is the tracer concentration at time t, location x, v is the drift velocity and D is related to the diffusivity of the media [4, 15]. The probabilistic approach to describe this flow from a mesoscopic view is given by the hypothesis that the path of a randomly chosen tracer particle is a Markov process that solves a stochastic differential equation (SDE) driven by Brownian motion. The basis of this hypothesis is that the conditional probability density function of this Markov process solves a forward equation of the same form as the ADE [7, 6].
However, for some heavy-tailed flows, the second order diffusion model can be inadequate. For such cases a model called the fractional advection-dispersion equation or fADE of the following form has been proposed [27]:
For the current discussion, we will consider a one-dimensional concentration $c(x,t)$ where x denotes the distance from the origin of the plume, v is the drift velocity and D is a function that changes with the diffusion of the tracer. The fractional differentiation order $\alpha \in (1,2)$ controls the tail of the flow. Here the fractional derivative of order α for any function f is defined as in [3] by:
(1.1)
\[ \frac{\partial c}{\partial t}=-\frac{\partial }{\partial x}\big[v(x)c\big]+\frac{{\partial ^{\alpha -1}}}{\partial {x^{\alpha -1}}}\bigg[D(x)\frac{\partial c}{\partial x}\bigg].\]
\[ \frac{{d^{\alpha }}f(x)}{d{x^{\alpha }}}=\frac{-1}{\varGamma (1-\alpha )}{\underset{0}{\overset{\infty }{\int }}}\big[f(x-u)-f(x)+u{f^{\prime }}(x)\big]\alpha {y^{-1-\alpha }}du.\]
The negative fractional derivative is given by $\frac{{d^{\alpha }}f(x)}{d{(-x)^{\alpha }}}=\frac{{d^{\alpha }}g(x)}{d{x^{\alpha }}}$, where $g(x)=f(-x)$.If we assume that the diffusion coefficient D is location invariant, then the fADE can be associated with an SDE driven by an α-stable Lévy process ${X_{t}}$ 1 of the form:
Following the Brownian diffusion argument, in the heavy-tailed diffusion model, we assume a random particle’s position at time t is given by the process ${Y_{t}}$ that solves SDE (1.2). It can be shown that ${Y_{t}}$ is a Markov process [1]. Let us denote the transition probability density function of ${Y_{t}}$ by ${p_{{y_{0}}}}(y,t)$, i.e. $P({Y_{t}}\in A|{Y_{0}}={y_{0}})={\int _{A}}{p_{{y_{0}}}}(y,t)dy$. For the initial distribution $\mu (u)$ the pdf of ${Y_{t}}$ is given by
e.g., for a ground water tracer concentration modeling it is reasonable to assume that all particles start at location ${y_{0}}$ at time $t=0$ hence $\mu (u)=I(u={y_{0}})$. If $1<\alpha <2$, it can be shown [10] that $p(y,t)$ solves the forward equation:
Typically the coefficient functions a and b depend on t via ${Y_{t}}$ as in (1.2) hence the terms in the forward equation are expressed as $a(y)$ and $b(y)$, where y is supposed to be in the range space of ${Y_{t}}$. However, without loss of generality, the functions can also be considered as a mapping from $T\times \mathbb{R}\to \mathbb{R}$ as $a(y,t)$ and $b(y,t)$ and the associated forward equation will remain the same.
(1.4)
\[\begin{aligned}{}\frac{\partial p(y,t)}{\partial t}=-\frac{\partial }{\partial y}\big[a(y)p(y,t)\big]& +\frac{(1+\beta )}{2}{\bigg[\hspace{-0.1667em}-\cos \bigg(\frac{\pi \alpha }{2}\bigg)\hspace{-0.1667em}\bigg]^{-1}}\hspace{-0.1667em}\hspace{-0.1667em}\cdot \frac{{\partial ^{\alpha }}}{\partial {y^{\alpha }}}\big[{b^{\alpha }}(y)p(y,t)\big]\\ {} & +\frac{(1-\beta )}{2}{\bigg[\hspace{-0.1667em}-\cos \bigg(\frac{\pi \alpha }{2}\bigg)\hspace{-0.1667em}\bigg]^{-1}}\hspace{-0.1667em}\hspace{-0.1667em}\cdot \frac{{\partial ^{\alpha }}}{\partial {(-y)^{\alpha }}}\big[{b^{\alpha }}(y)p(y,t)\big].\end{aligned}\]Thus for a location invariant D, the fADE in (1.1) can be written as:
which essentially has the same form as the forward equation given in (1.4) with $\beta =1$, $a(x,t)=v(x,t)$, ${[-\cos (\frac{\pi \alpha }{2})]^{-1}}{b^{\alpha }}(x,t)\equiv D(t)$. This shows that in a heavy-tailed plume that follows fADE (1.5), the position ${Y_{t}}$ of a randomly chosen tracer particle at time t solves SDE (1.2). Choosing $\beta \ne 1$ will allow the particle transition in forward or backward direction.
(1.5)
\[ \frac{\partial C}{\partial t}=-\frac{\partial }{\partial x}\big[v(x,t)C\big]+D(t)\frac{{\partial ^{\alpha }}}{\partial {x^{\alpha }}}\bigg[\frac{\partial C}{\partial x}\bigg],\]In most practical tracking scenarios, the plume data consists of observed tracer concentration values $c(x,t)$ over a range of locations at certain time points. We write $p(y,t)={K_{t}}c(y,t)$, where ${K_{t}}$ is a suitable scale parameter that adjusts the total mass for $c(x,t)$ so that ${\int _{\mathbb{R}}}p(y,t)dy=1$ for each t.
In case $v(\cdot )$ and $D(\cdot )$ are constants, (1.2) reduces to ${Y_{t}}=at+b{X_{t}}$ with $a=v$, $b={[-\cos (\frac{\pi \alpha }{2})D]^{1/\alpha }}$ and consequently, ${Y_{t}}$ itself is a stable process. The parameters of this stable process then can be estimated using the observed data (see [11] for details). However, when the coefficient functions are not constants, the solution to the SDE in (1.2) may not have a closed form.
In this paper, we present a numerical method to solve ${Y_{t}}$ through the fitted fADE equation using the observed tracer concentration data. Also, often for a groundwater contamination modeling problem the functional form of the advection or dispersion coefficients are not known to start with. In that case, the proper advection and dispersion coefficient functions are needed to be estimated. This presents an inverse problem. We formulate this inverse problem as an optimization problem and develop a Levenberg–Marquardt (L-M) regularization method to obtain the proper advection and dispersion coefficient function from the observed data (the case for a linear drift function is presented here). Then we proceed to generate ${Y_{t}}$ using these fitted coefficient functions. The detailed methodology is described in Section 2. Section 3 provides an illustrative example. Section 4 includes concluding remarks.
2 Methodology
We propose a method for generating a weak solution for the SDE given in (1.2). The core technique is the standard probability integral transformation. Suppose the probability density function (pdf) of the solution process ${Y_{t}}$ at a given time point t is given by $p(y,t)$ and the cumulative distribution function (CDF) is defined by ${F_{{Y_{t}}}}(x)={\int _{-\infty }^{x}}p(y,t)dy$. Since $F({Y_{t}})\sim \mathit{Uniform}(0,1)$ [21], a random observation of ${Y_{t}}$, say y, can be generated by ${F_{{Y_{t}}}^{-1}}(u)=y$ where u is a random observation from $\mathit{Uniform}(0,1)$. For the present problem, closed form of $p(t,y)$ is not available, and we use point-wise numerical approximations to estimate $p(y,t)$ for any fixed t and fixed y. Therefore we can only use a numerical scheme to approximate ${F_{{Y_{t}}}}$ and ${F_{{Y_{t}}}^{-1}}$.
2.1 Simulating a solution process for an SDE driven by a heavy-tailed stable process
For a heavy-tailed plume where the particle path can be modeled by the SDE in (1.2), the forward equation is of the form (1.5). To generate the solution to this SDE we can consider two scenarios: (i) we may want to solve an SDE with a given form of coefficient functions $v(\cdot )$ and $D(\cdot )$ and driven by a given α-stable process; (ii) a more practical problem when only the tracer concentrations at different time points over a range of locations are observed from a plume, while the coefficient functions need to be estimated along with the parameters of the driving process ${X_{t}}$ in (1.2) from the observed concentration values.
In either case, first, we find $p(t,y)$ that solves the forward equation (1.4) with given or estimated coefficient functions and then generate ${Y_{t}}$ using these $p(t,y)$’s.
${Y_{t}}$ Simulation Steps:
Step 1: Numerical Solution for $p(t,y)$:
Assuming the probability density $p(y,t)={K_{t}}c(y,t)$, we estimate ${K_{t}}$, α and β from the observed data that assumed to satisfy fADE (1.5). $p(y,t)$ can be solved numerically using a discretization scheme from the forward equation (1.4). In case the coefficient functions in (1.4) are not known, an inverse problem optimization technique can be used for estimating the parameters. See Section 2.2 for details of this step.
In either case the discretization method provides only point-wise approximation of $p(y,t)$ values for given t and y and not a closed functional form.
Step 2: Generating ${Y_{t}}$:
To generate ${Y_{t}}$ for a fixed t using an inverse CDF transformation, we present a straightforward scheme with a trapezoidal rule to approximate ${F_{{Y_{t}}}}$ and ${F_{{Y_{t}}}^{-1}}$. Approximate ${Y_{t}}$’s can be simulated using $p(y,t)$’s from step 1 as follows:
-
(a) Choose suitable grid ${y_{0}}<{y_{1}}<\cdots <{y_{n}}$ and solve for $p({y_{i}},t)$ as mentioned in step 1. The endpoints ${y_{0}}$ and ${y_{n}}$ can be chosen so that $p(y,t)$ is small or negligible for $y<{y_{0}}$ or $y>{y_{n}}$. The fitted coefficient and parameters from the observed data in step 1 can be used for any ${y_{i}}$’s.
-
(b) Approximate ${F_{{Y_{t}}}}({y_{i}})$ by ${\hat{F}_{{Y_{t}}}}({y_{i}})$ for $i=1,2,\dots ,n$, using the trapezoidal rule:\[\begin{aligned}{}{\hat{F}_{{Y_{t}}}}({y_{1}})& =({y_{1}}-{y_{0}})\frac{p({y_{1}},t)+p({y_{0}},t)}{2},\\ {} {\hat{F}_{{Y_{t}}}}({y_{2}})& ={\hat{F}_{{Y_{t}}}}({y_{1}})+({y_{2}}-{y_{1}})\frac{p({y_{2}},t)+p({y_{1}},t)}{2},\\ {} \vdots \\ {} {\hat{F}_{{Y_{t}}}}({y_{r}})& ={\hat{F}_{{Y_{t}}}}({y_{r-1}})+({y_{r}}-{y_{r-1}})\frac{p({y_{r}},t)+p({y_{r-1}},t)}{2},\\ {} \vdots \\ {} {\hat{F}_{{Y_{t}}}}({y_{n}})& ={\hat{F}_{{Y_{t}}}}({y_{n-1}})+({y_{n}}-{y_{n-1}})\frac{p({y_{n}},t)+p({y_{n-1}},t)}{2}.\end{aligned}\]
-
(c) For a randomly generated $\mathit{Uniform}(0,1)$ observation u, find its placement in the grid, i.e. ${y_{r}}$, ${y_{r+1}}$ such that ${\hat{F}_{{Y_{t}}}}({y_{r}})\le u<{\hat{F}_{{Y_{t}}}}({y_{r+1}})$, for some $r\in \{0,1,\dots ,n\}$. A way of generating associated random observation of ${Y_{t}}$ or ${F_{{Y_{t}}}^{-1}}(u)$ (approximated) as $\hat{y}$ is described below:\[ \hat{y}=\left\{\begin{array}{l@{\hskip10.0pt}l}{y_{r}},& if\hspace{2.5pt}u={y_{r}};\\ {} {y_{r}}+\frac{u-{\hat{F}_{{Y_{t}}}}({y_{r}})}{p(t,{y_{r}})},& if\hspace{2.5pt}{y_{r}}<u<{y_{r+1}}.\end{array}\right.\]In case $u<{\hat{F}_{{Y_{t}}}}({y_{0}})$ set $\hat{y}={y_{0}}$ and if $u\ge {\hat{F}_{{Y_{t}}}}({y_{n}})$ set $y={y_{n}}$.
-
(d) To generate N i.i.d approximated random observations ${\hat{y}_{1}},{\hat{y}_{2}},\dots ,{\hat{y}_{N}}$ of ${Y_{t}}$, start with randomly generated independent $\mathit{Uniform}(0,1)$ observations ${u_{1}},{u_{2}},\dots {u_{N}}$ and then repeat step (c). Note that the efficiency of this approximation will depend on the spacing of the chosen grid since we are actually generating sample from distribution ${\hat{F}_{{Y_{t}}}}$ here and ${\hat{F}_{{Y_{t}}}}\to {F_{{Y_{t}}}}$ as $\underset{0\le i\le n}{\max }|{y_{i+1}}-{y_{i}}|\to 0$ following the convergence property of the standard trapezoidal rule.
2.2 A formulation for a numerical solution of the forward equation
In this subsection, we will present a finite volume scheme with the characteristic line for solving a forward equation of the form (1.4). To incorporate a discretized finite volume scheme with proper boundary condition, we re-parameterize the forward equation with fractional order $2-\lambda $ with $0<\lambda <1$ (consistent with the construction described in [25]) as follows:
Essentially, compared to (1.4) the fractional derivative order $\alpha =2-\lambda $, while $\gamma =\frac{1+\beta }{2}$ indicates the relative weight of forward versus backward transition. ${_{{x_{l}}}}{D_{x}^{-\lambda }}$ and ${_{x}}{D_{{x_{r}}}^{-\lambda }}$ represent the left and right fractional integral operators defined as
To comply with the boundary conditions, $x={x_{l}}$ and $x={x_{r}}$ are set to be the inflow and outflow boundaries respectively [25], while $a({x_{l}})$ and $a({x_{r}})$ are assumed to be non-negative.
(2.1)
\[ \begin{aligned}{}& \frac{\partial p}{\partial t}+\frac{\partial }{\partial x}a(x)p-\frac{\partial }{\partial x}\big[b(t)\big({\gamma _{{x_{l}}}}{D_{x}^{-\lambda }}+{(1-\gamma )_{x}}{D_{{x_{r}}}^{-\lambda }}\big)Dp\big]=0,\hspace{1em}{x_{l}}<x<{x_{r}},\\ {} & p({x_{l}},t)=0,\hspace{2em}p({x_{r}},t)=0,\hspace{1em}0<t<T.\end{aligned}\](2.2)
\[ \begin{aligned}{}{_{{x_{l}}}}{D_{x}^{-\lambda }}f(x):& =\frac{1}{\varGamma (\beta )}{\int _{{x_{l}}}^{x}}{(x-s)^{\lambda -1}}f(s)\mathrm{d}s,\\ {} {_{x}}{D_{{x_{r}}}^{-\lambda }}f(x):& =\frac{1}{\varGamma (\beta )}{\int _{x}^{{x_{r}}}}{(s-x)^{\lambda -1}}f(s)\mathrm{d}s.\end{aligned}\]2.2.1 The discretized finite volume scheme and the accumulation terms
Let us define a partition ${\mathcal{P}_{t}}$ on the time interval $[0,T]$, and a partition ${\mathcal{P}_{x}}$ on the space interval $[{x_{l}},{x_{r}}]$
with $\Delta {t_{n}}={t_{n}}-{t_{n-1}}$ for $n=1,\dots ,N$ and $\Delta {x_{i}}={x_{i}}-{x_{i-1}}$ for $i=1,\dots ,I$. Let $\bar{{\mathcal{P}_{x}}}$ be the dual partition of ${\mathcal{P}_{x}}$ defined by ${x_{i-1/2}}=\frac{1}{2}({x_{i-1}}+{x_{i}})$ with $\Delta {x_{i+1/2}}={x_{i+1/2}}-{x_{i-1/2}}$ for $i=1,\dots ,I-1$ as well as $\Delta {x_{1/2}}={x_{1/2}}-{x_{0}}$ and $\Delta {x_{I+1/2}}={x_{I}}-{x_{I-1/2}}$.
(2.3)
\[ \begin{aligned}{}{\mathcal{P}_{t}}:0& ={t_{0}}<{t_{1}}<\cdots <{t_{N}}=T,\\ {} {\mathcal{P}_{x}}:{x_{l}}& ={x_{0}}<{x_{1}}<\cdots <{x_{I}}={x_{r}},\end{aligned}\]Let $y=r(\theta ;\bar{x},\bar{t})$ be a continuous and piecewise-smooth curve that passes through the point $\bar{x}$ at time $\bar{t}$, such that $r(\theta ;{x_{i-1/2}},{t_{n}})$ and $r(\theta ;{x_{i+1/2}},{t_{n}})$ do not meet each other during the time period $[{t_{n-1}},{t_{n}}]$. We define a space-time control volume ${\varOmega _{i}^{n}}$ by extending the cell $[{x_{i-1/2}},{x_{i+1/2}}]$ along the curve $r(t;x,{t_{n}})$ from $t={t_{n}}$ to $t={t_{n-1}}$
Assuming the prism ${\varOmega _{i}^{n}}$ does not intersect the boundary $x={x_{l}}$ and $x={x_{r}}$ of the domain during the time period $[{t_{n-1}},{t_{n}}]$, let ${x^{\ast }}=r({t_{n-1}};x,{t_{n}})$ be the foot at time ${t_{n-1}}$ of the curve with head x at time ${t_{n}}$.
(2.4)
\[ {\varOmega _{i}^{n}}=\big\{(x,t):r(t;{x_{i-1/2}},{t_{n}})<x<r(t;{x_{i+1/2}},{t_{n}}),\hspace{1em}{t_{n-1}}<t<{t_{n}}\big\}.\]Integrating equation (2.1) over the control volume ${\varOmega _{i}^{n}}$ we get
Without loss of generality, the accumulation term in (2.4) can be evaluated by assuming
and accumulation term can be re-written as follows:
Here the notation $t(x;{x_{i-1/2}},{t_{n}})$ represents the time instant that the curve $r(t;,{x_{i-1/2}},{t_{n}})=x$. The notation $t(x;{x_{i+1/2}},{t_{n}})$ is defined similarly.
(2.5)
\[ {\int _{{\varOmega _{i}^{n}}}}\frac{\partial p}{\partial t}dxdt+{\int _{{t_{n-1}}}^{{t_{n}}}}\Big(a(x)p-b(t)\frac{{\partial ^{1-\lambda }}}{\partial {x^{1-\lambda }}}p\Big){\Big|_{r(t;{x_{i-1/2}},{t_{n}})}^{r(t;{x_{i+1/2}},{t_{n}})}}dt=0.\](2.7)
\[\begin{aligned}{}{\int _{{\varOmega _{i}^{n}}}}\frac{\partial p}{\partial t}dxdt& ={\int _{{x_{i-1/2}^{\ast }}}^{{x_{i-1/2}}}}\Bigg[{\int _{{t_{n-1}}}^{t(x;{x_{i-1/2}},{t_{n}})}}\frac{\partial p}{\partial t}\hspace{0.2222em}dt\Bigg]\hspace{0.2222em}dx\\ {} & \hspace{1em}+{\int _{{x_{i-1/2}}}^{{x_{i+1/2}^{\ast }}}}\Bigg[{\int _{{t_{n-1}}}^{{t_{n}}}}\frac{\partial p}{\partial t}\hspace{0.2222em}dt\Bigg]\hspace{0.2222em}dx+{\int _{{x_{i+1/2}^{\ast }}}^{{x_{i+1/2}}}}\Bigg[{\int _{t(x;{x_{i+1/2}},{t_{n}})}^{{t_{n}}}}\frac{\partial p}{\partial t}\hspace{0.2222em}dt\Bigg]\hspace{0.2222em}dx.\end{aligned}\]A simple calculation of the second term on the right-hand side of (2.7) yields
The first and third terms on the right-hand side of (2.7) are integrated as follows
and
(2.8)
\[ {\int _{{x_{i-1/2}}}^{{x_{i+1/2}^{\ast }}}}\Bigg[{\int _{{t_{n-1}}}^{{t_{n}}}}\frac{\partial p}{\partial t}\hspace{0.2222em}dt\Bigg]dx={\int _{{x_{i-1/2}}}^{{x_{i+1/2}^{\ast }}}}p(x,{t_{n}})dx-{\int _{{x_{i-1/2}}}^{{x_{i+1/2}^{\ast }}}}p(x,{t_{n-1}})dx.\](2.9)
\[\begin{aligned}{}& {\int _{{x_{i-1/2}^{\ast }}}^{{x_{i-1/2}}}}\Bigg[{\int _{{t_{n-1}}}^{t(x;{x_{i-1/2}},{t_{n}})}}\frac{\partial p}{\partial t}\hspace{0.2222em}dt\Bigg]dx\\ {} & \hspace{1em}={\int _{{t_{n-1}}}^{{t_{n}}}}p\big(r(t;{x_{i-1/2}},{t_{n}}),t\big)\frac{\partial r(t;{x_{i-1/2}},{t_{n}})}{\partial t}\hspace{0.2222em}dt-{\int _{{x_{i-1/2}^{\ast }}}^{{x_{i-1/2}}}}p(x,{t_{n-1}})\hspace{0.2222em}dx,\end{aligned}\](2.10)
\[\begin{aligned}{}& {\int _{{x_{i+1/2}^{\ast }}}^{{x_{i+1/2}}}}\Bigg[{\int _{t(x;{x_{i+1/2}},{t_{n}})}^{{t_{n}}}}\frac{\partial p}{\partial t}\hspace{0.2222em}dx\Bigg]dt\\ {} & \hspace{1em}={\int _{{x_{i+1/2}^{\ast }}}^{{x_{i+1/2}}}}p(x,{t_{n}})\hspace{0.2222em}dx-{\int _{{t_{n-1}}}^{{t_{n}}}}p\big(r(t;{x_{i+1/2}},{t_{n}}),t\big)\frac{\partial r(t;{x_{i+1/2}},{t_{n}})}{\partial t}\hspace{2.5pt}dt.\end{aligned}\]Incorporating equations (2.8)–(2.18) into (2.7) we get
The derivation shows that (2.11) does not depend on the assumption (2.4).
(2.11)
\[ \begin{aligned}{}{\int _{{\varOmega _{i}^{n}}}}\frac{\partial p}{\partial t}dxdt& ={\int _{{x_{i-1/2}}}^{{x_{i+1/2}}}}p(x,{t_{n}})dx-{\int _{{x_{i-1/2}^{\ast }}}^{{x_{i+1/2}^{\ast }}}}p(x,{t_{n-1}})dx\\ {} & \hspace{1em}+{\int _{{t_{n-1}}}^{{t_{n}}}}p\big(r(t;{x_{i-1/2}},{t_{n}}),t\big)\frac{\partial r(t;{x_{i-1/2}},{t_{n}})}{\partial t}\hspace{0.2222em}dt\\ {} & \hspace{1em}-{\int _{{t_{n-1}}}^{{t_{n}}}}p\big(r(t;{x_{i+1/2}},{t_{n}}),t\big)\frac{\partial r(t;{x_{i+1/2}},{t_{n}})}{\partial t}\hspace{0.2222em}dt.\end{aligned}\]Now we just set the space-time boundaries $r(t;{x_{i\pm 1/2}},{t_{n}})$ of the control volume ${\varOmega _{i}^{n}}$ to be the characteristic curves, which are defined by the initial value problem of the ordinary differential equation
That is, because the characteristics $r(t;{x_{i\pm 1/2}},{t_{n}})$ are assumed to be tracked exactly and hence the residual advection term vanishes naturally.
Substituting (2.11) for the accumulation term in (2.5), we obtain a locally conservative reference equation on an interior space-time control volume ${\varOmega _{i}^{n}}$ as follows:
which can be approximated as
with the Lagrangian interface fluxes
(2.13)
\[ \begin{aligned}{}& {\int _{{x_{i-1/2}}}^{{x_{i+1/2}}}}p(x,{t_{n}})\hspace{0.2222em}dx+{\int _{{t_{n-1}}}^{{t_{n}}}}F(p,b)\big(r(t;{x_{i+1/2}},{t_{n}}),t\big)\hspace{0.2222em}dt\\ {} & \hspace{2em}-{\int _{{t_{n-1}}}^{{t_{n}}}}F(p,b)\big(r(t;{x_{i-1/2}},{t_{n}}),t\big)\hspace{0.2222em}dt\\ {} & \hspace{1em}={\int _{{x_{i-1/2}^{\ast }}}^{{x_{i+1/2}^{\ast }}}}p(x,{t_{n-1}})\hspace{2.5pt}dx,\end{aligned}\]Note that the interface fluxes $F(p,b)(r(t;{x_{i\pm 1/2}},{t_{n}}),t)$ are defined across the space-time boundary $r(t;{x_{i-1/2}},{t_{n}})$ and $r(t;{x_{i+1/2}},{t_{n}})$ of the space-time control volume ${\varOmega _{i}^{n}}$.
To discretize equation (2.1), we further let ${\{{\phi _{i}}\}_{i=1}^{I}}$ be a set of hat functions such that ${\phi _{i}}({x_{i}})=1$ and ${\phi _{i}}({x_{j}})=0$ for $j\ne i$. Hence the finite volume approximation ${p_{h}}$ to the true solution p can be expressed as
and the finite volume scheme on the interval $[{x_{i-1/2}},{x_{i+1/2}}]$ for $i=1,\dots ,I-1$ has the form
where, ${[{z_{i,j}}]_{i,j=1}^{I-1}}$ is the coefficient matrix of the fractional term, and is given by
(2.16)
\[ {\int _{{x_{i-1/2}}}^{{x_{i+1/2}}}}p(x,{t_{n}})\hspace{0.2222em}dx+\Delta t{\sum \limits_{j=1}^{I-1}}{p_{j}^{n}}{z_{i,j}}={\int _{{x_{i-1/2}^{\ast }}}^{{x_{i+1/2}^{\ast }}}}p(x,{t_{n-1}})\hspace{2.5pt}dx,\](2.17)
\[ \begin{array}{r}\displaystyle {z_{i,j}}=b\big[\big({\gamma _{0}}{D_{{x_{i-1/2}}}^{-\lambda }}+{(1-\gamma )_{{x_{1-1/2}}}}{D_{1}^{-\lambda }}\big)D{\phi _{j}}({x_{i-1/2,{t_{n}}}})\\ {} \displaystyle \hspace{2em}\hspace{1em}-\big({\gamma _{0}}{D_{{x_{+1/2}}}^{-\lambda }}+{(1-\gamma )_{{x_{1+1/2}}}}{D_{1}^{-\lambda }}\big)D{\phi _{j}}({x_{i+1/2,{t_{n}}}})\big].\end{array}\]Assuming the trial functions $p(x,{t_{n}})$ are chosen to be piecewise linear functions on $[{x_{l}},{x_{r}}]$ with respect to the fixed spatial partition ${\mathcal{P}_{x}}$, we evaluate the accumulation term at time step ${t_{n}}$ in the finite volume scheme (2.16) analytically as
In addition, we can compute ${x_{i-1/2}^{\ast }}$ and ${x_{i+1/2}^{\ast }}$ at time step ${t_{n-1}}$ via a backward approximate characteristic tracking. Since the trial function $p(x,{t_{n-1}})$ is also piecewise linear with respect to the fixed spatial partition ${\mathcal{P}_{x}}$ at time ${t_{n-1}}$, we can evaluate the accumulation term at time step ${t_{n-1}}$ analytically. Because the accumulation term at time step ${t_{n-1}}$ affects only the right-hand side of the finite volume scheme (2.16), the scheme retains a symmetric and positive-definite coefficient matrix. Furthermore, the finite volume scheme (2.16) is locally conservative, even if the characteristics are computed approximately.
2.2.2 The fractional diffusion term and the stiffness matrix
The finite volume scheme (2.16) appears similar to the one for the canonical second-order diffusion equation, but has a fundamental difference. Although the hat functions ${\phi _{j}}$ have local support, (2.2) reveals that ${_{{x_{l}}}}{D_{x}^{-\lambda }}{\phi _{j}}$ and ${_{x}}{D_{{x_{r}}}^{-\lambda }}{\phi _{j}}$ have global support. Therefore, the stiffness matrix $Z:={[{z_{i,j}}]_{i,j=1}^{I-1}}$ is a full matrix, which requires $O({N^{2}})$ of memory to store. Numerical schemes for space-fractional differential equations were traditionally solved by the Gaussian type direct solvers that require $O({N^{3}})$ of computations [13, 14, 18, 19]. In recent years, there have been some other notable developments in methods for solving the algebraic linear systems arising from discretization of fractional-order problems, especially for space dimension higher than one (see [9]); we plan to explore them in our future work.
To simplify the computation of the diffusion term, we have to explore the structure of the stiffness matrix Z [17, 26]. The entries of the stiffness matrix Z are presented below.
Its diagonal entries are given by:
for $1\le i\le I-1$.
(2.19)
\[ \begin{aligned}{}{z_{i,i}}& =\frac{1}{\varGamma (\lambda +1){2^{\lambda }}}b\gamma {h^{\lambda -1}}+\frac{1}{\varGamma (\lambda +1)}b(1-\gamma ){h^{\lambda -1}}\bigg(2{\left(\frac{1}{2}\right)^{\lambda }}-{\left(\frac{3}{2}\right)^{\lambda }}\bigg)\\ {} & \hspace{1em}+\frac{1}{\varGamma (\lambda +1)}b\gamma {h^{\lambda -1}}\bigg(2{\left(\frac{1}{2}\right)^{\lambda }}-{\left(\frac{3}{2}\right)^{\lambda }}\bigg)\\ {} & \hspace{1em}+\frac{1}{\varGamma (\lambda +1){2^{\lambda }}}b(1-\gamma ){h^{\lambda -1}}\\ {} & =\frac{1}{\varGamma (\lambda +1)}b{h^{\lambda -1}}\bigg({2^{-\lambda }}+2{\left(\frac{1}{2}\right)^{\lambda }}-{\left(\frac{3}{2}\right)^{\lambda }}\bigg),\end{aligned}\]Whereas, the sub-triangular entries of the matrix Z are given by:
for $2\le i\le I-1$, and
for $3\le i\le I-1$ and $1\le j\le i-2$.
(2.20)
\[ \begin{aligned}{}{z_{i,i-1}}& =-\frac{1}{\varGamma (\lambda +1)}b\gamma {h^{\lambda -1}}\bigg(2{\left(\frac{1}{2}\right)^{\lambda }}-{\left(\frac{3}{2}\right)^{\lambda }}\bigg)-\frac{1}{\varGamma (\lambda +1){2^{\lambda }}}b(1-\gamma ){h^{\lambda -1}}\\ {} & \hspace{1em}-\frac{1}{\varGamma (\lambda +1)}b\gamma {h^{\lambda -1}}\bigg({\left(\frac{5}{2}\right)^{\lambda }}-2{\left(\frac{3}{2}\right)^{\lambda }}+{\left(\frac{1}{2}\right)^{\lambda }}\bigg)\\ {} & =\frac{1}{\varGamma (\lambda +1)}b{h^{\lambda -1}}\bigg(3{\left(\frac{3}{2}\right)^{\lambda }}\gamma -3{\left(\frac{1}{2}\right)^{\lambda }}\gamma -{\left(\frac{5}{2}\right)^{\lambda }}\gamma -{2^{-\lambda }}(1-\gamma )\bigg),\end{aligned}\](2.21)
\[ \begin{aligned}{}{z_{i,j}}& =\frac{1}{\varGamma (\lambda +1)}b\gamma {h^{\lambda -1}}\bigg({\left(i-j+\frac{1}{2}\right)^{\lambda }}-2{\left(i-j-\frac{1}{2}\right)^{\lambda }}+{\left(i-j-\frac{3}{2}\right)^{\lambda }}\bigg)\\ {} & \hspace{1em}-\frac{1}{\varGamma (\lambda +1)}b\gamma {h^{\lambda -1}}\bigg({\left(i-j+\frac{3}{2}\right)^{\lambda }}-2{\left(i-j+\frac{1}{2}\right)^{\lambda }}+{\left(i-j-\frac{1}{2}\right)^{\lambda }}\bigg)\end{aligned}\]The super-triangular entries of matrix Z can be also derived as
for $1\le i\le I-2$, and
for $1\le j\le I-3$ and $j+2\le i\le I-1$.
(2.22)
\[\begin{aligned}{}{z_{i,i+1}}& =-\frac{1}{\varGamma (\lambda +1)}b(1-\gamma ){h^{\lambda -1}}\bigg({\left(\frac{1}{2}\right)^{\lambda }}-2{\left(\frac{3}{2}\right)^{\lambda }}+{\left(\frac{5}{2}\right)^{\lambda }}\bigg)\\ {} & \hspace{1em}-\frac{1}{\varGamma (\lambda +1){2^{\lambda }}}b\gamma {h^{\lambda -1}}-\frac{1}{\varGamma (\lambda +1)}b(1-\gamma ){h^{\lambda -1}}\bigg(2{\left(\frac{1}{2}\right)^{\lambda }}-{\left(\frac{3}{2}\right)^{\lambda }}\bigg)\\ {} & =\frac{1}{\varGamma (\lambda +1)}b{h^{\lambda -1}}\bigg(3{\left(\frac{3}{2}\right)^{\lambda }}(1-\gamma )-3{\left(\frac{1}{2}\right)^{\lambda }}(1-\gamma )\\ {} & \hspace{1em}-{\left(\frac{5}{2}\right)^{\lambda }}(1-\gamma )-{2^{-\lambda }}\gamma \bigg),\end{aligned}\](2.23)
\[\begin{aligned}{}{z_{i,j}}& =\frac{1}{\varGamma (\lambda +1)}b(1-\gamma ){h^{\lambda -1}}\bigg(2{\left(j-i+\frac{1}{2}\right)^{\lambda }}\\ {} & \hspace{1em}-{\left(j-i-\frac{1}{2}\right)^{\lambda }}-{\left(j-i+\frac{3}{2}\right)^{\lambda }}\bigg)\\ {} & \hspace{1em}-\frac{1}{\varGamma (\lambda +1)}b(1-\gamma ){h^{\lambda -1}}\bigg(2{\left(j-i-\frac{1}{2}\right)^{\lambda }}\\ {} & \hspace{1em}-{\left(j-i-\frac{3}{2}\right)^{\lambda }}-{\left(j-i+\frac{1}{2}\right)^{\lambda }}\bigg)\end{aligned}\]2.2.3 Estimation of coefficient function
We present the case where the drift coefficient function in fADE (1.5) (and therefore in equation (2.1)) is assumed to be a piecewise linear function of x,
where parameters ${a_{0}}$, ${a_{1}}$, ${a_{2}}$ and ${a_{3}}$ are to be estimated from the observed concentration data. The main idea on this part is to obtain certain measurements through physical or mechanical experiments, and then use the data to calibrate these parameters in the fADE (see [16]). This is an inverse problem: based on the initial guess ${p_{0}}$ of the equation (2.1), and certain observation (or desired) data such as values of the state variable g at the final time, we attempt to seek for the constant parameters ${a_{0}}$, ${a_{1}}$, ${a_{2}}$ and ${a_{3}}$ from the governing differential equation (2.1).
(2.24)
\[ a(x)=\left\{\begin{array}{l@{\hskip10.0pt}l}{a_{0}}-{a_{1}}x,& \text{if}\hspace{2.5pt}{x_{l}}\le x\le {x_{m}},\\ {} {a_{2}}-{a_{3}}x,& \text{if}\hspace{2.5pt}{x_{m}}<x\le {x_{r}},\end{array}\right.\]We formulate the inverse problem as an optimization and develop a Levenberg–Marquardt (L-M) regularization method (see, [12, 20, 24]) to iteratively identify the parameter. It is known that the inverse problem usually requires multiple runs of the forward problem. Considering the computational cost of the forward problem is already high, the inverse problem could become infeasible. Hence we propose an optimization algorithm for the parameter estimation.
Here we only present the details of the fitting of ${a_{0}}$ and ${a_{1}}$. The other two parameters ${a_{2}}$ and ${a_{3}}$ can be estimated similarly. The parameter identification of ${a_{0}}$ and ${a_{1}}$ can be formulated as follows: let $\boldsymbol{\alpha }:=\{{a_{0}},{a_{1}}\}$, then to find ${\alpha _{inv}}$ that satisfies
where, $\hat{p}({x_{i}},t)={\hat{K}_{t}}c({x_{i}},t)$ for observed concentration c at location ${x_{i}}$, time t.
(2.25)
\[ {\alpha _{inv}}=\text{arg}\hspace{2.5pt}\hspace{2.5pt}\underset{\boldsymbol{\alpha }}{\min }\hspace{2.5pt}\hspace{2.5pt}\mathcal{G}(\boldsymbol{\alpha }):=\frac{1}{2}{\sum \limits_{i=1}^{M}}{\big[p({x_{i}},t;\boldsymbol{\alpha })-\hat{p}({x_{i}},t)\big]^{2}},\]In case the data is available for time points ${t_{1}},{t_{2}},\dots ,{t_{R}}$ we rewrite equation (2.25) as
here ${x_{{t_{k}},i}}$ is the ith observed location at time ${t_{k}}$ and ${w_{k}}$ is the weight assigned for the sample set available for the same time point. (For example the data set used in Section 3 has $R=3$).
(2.26)
\[ \begin{aligned}{}{\alpha _{inv}}& =\text{arg}\hspace{2.5pt}\hspace{2.5pt}\underset{\boldsymbol{\alpha }}{\min }\hspace{2.5pt}\hspace{2.5pt}\mathcal{G}(\boldsymbol{\alpha })\\ {} & =\frac{1}{2}{\sum \limits_{k=1}^{R}}{\sum \limits_{i=1}^{{M_{k}}}}{w_{k}}{\big[p({x_{{t_{k}},i}},{t_{k}};\boldsymbol{\alpha })-\hat{p}({x_{{t_{k}}i}},{t_{k}})\big]^{2}},\end{aligned}\]An iteration algorithm such as the Newton method with line searching could be employed to find the solution to (2.26). Basically, the Newton algorithm for minimizing (2.26) uses the first and second derivatives of the objective function $\mathcal{G}(\boldsymbol{\alpha })$:
where k represents the kth iteration. It is easy to check that (2.27) is equivalent to solve
where
and
for $j=1,\dots ,{M_{i}}$, and $i=1,2,\dots ,R$. Note that in practice, we always use the finite difference
(2.27)
\[ {\boldsymbol{\alpha }_{k+1}}={\boldsymbol{\alpha }_{k}}-\frac{{\mathcal{G}^{\prime }}({\boldsymbol{\alpha }_{k}})}{{\mathcal{G}^{\prime\prime }}({\boldsymbol{\alpha }_{k}})},\](2.28)
\[ {\boldsymbol{\alpha }_{k+1}}={\boldsymbol{\alpha }_{k}}-{\big({\mathbf{J}_{k}^{T}}{\mathbf{J}_{k}}\big)^{-1}}{\mathbf{J}_{k}^{T}}{\mathbf{r}_{k}},\](2.29)
\[ \begin{aligned}{}{\mathbf{J}_{k}}& ={\big({\mathbf{J}^{1}};{\mathbf{J}^{2}};\dots ;{\mathbf{J}^{R}}\big)^{T}},\hspace{1em}\text{with}\\ {} {\mathbf{J}^{i}}& =\left(\frac{\partial p({x_{{t_{i}},1}},{t_{i}};\boldsymbol{\alpha })}{\partial \boldsymbol{\alpha }},\dots ,\frac{\partial p({x_{{t_{i}},{M_{i}}}},{t_{i}};\boldsymbol{\alpha })}{\partial \boldsymbol{\alpha }}\right),\end{aligned}\](2.30)
\[ \begin{aligned}{}{\mathbf{r}_{k}}& ={\big({\mathbf{r}^{1}};{\mathbf{r}^{2}};\dots ;{\mathbf{r}^{R}}\big)^{T}},\hspace{1em}\text{with}\\ {} {\mathbf{r}^{i}}& =({r_{i,1}},\dots ,{r_{i,{M_{i}}}}),\hspace{1em}{r_{i,j}}=p({x_{{t_{i}},j}},{t_{i}};\boldsymbol{\alpha })-\hat{p}({x_{{t_{i}},j}},{t_{i}}),\end{aligned}\]
\[ \frac{p({x_{{t_{i}},j}},{t_{i}};\boldsymbol{\alpha }+\delta )-p({x_{{t_{i}},j}},{t_{i}};\boldsymbol{\alpha })}{\delta }\]
with a small enough δ to approximate the derivatives in (2.29).However, the Newton method may fail to work because of ${\mathbf{J}_{k}^{T}}{\mathbf{J}_{k}}$ may be nearly zero. Therefore, the search direction ${d_{k}}:=-{\mathbf{J}_{k}^{T}}{\mathbf{r}_{k}}/{\mathbf{J}_{k}^{T}}{\mathbf{J}_{k}}$ may not point in a decent direction.
A common technique to overcome this kind of problem is the L-M algorithm (or Levenberg algorithm since a single parameter case is considered in this paper), which modifies (2.27) by the following formulation
where ${\varrho _{k}}$ is a positive penalty parameter, and ${I_{2}}$ is a $2\times 2$ identity matrix. The method coincides with the Newton algorithm when ${\varrho _{k}}=0$; and it gives a step closer to the gradient descent direction when ${\varrho _{k}}$ is large.
(2.31)
\[ {\boldsymbol{\alpha }_{k+1}}={\boldsymbol{\alpha }_{k}}-{\big({\mathbf{J}_{k}^{T}}{\mathbf{J}_{k}}+{\varrho _{k}}{I_{2}}\big)^{-1}}{\mathbf{J}_{k}^{T}}{\mathbf{r}_{k}},\]
Algorithm A: Parameter Identification Algorithm:
Given the observation data g and an initial guess ${\boldsymbol{\alpha }_{0}}={\{{a_{0}},{a_{1}}\}}_{0}$, choose $\rho \in (0,1)$, $\sigma \in (0,1/2)$, ${\varrho _{0}}>0$, and δ small enough. For $k=0,1,\dots ,M$, (let $M={M_{1}}+{M_{2}}+{M_{3}}$),
-
(1) Solve the equation (2.1) corresponding to ${\boldsymbol{\alpha }_{k}}$ and ${\boldsymbol{\alpha }_{k}}+\delta $ respectively to obtain $p(\cdot ,\cdot ,{\boldsymbol{\alpha }_{k}})$ and $p(\cdot ,\cdot ,{\boldsymbol{\alpha }_{k}}+\delta )$.
-
(3) Determine the search step ${\rho ^{m}}$ by Armijo rule: where m is the smallest non-negative integer.
-
(4) If $|{\rho ^{m}}{\mathbf{d}_{k}}|\le \text{Tolerance level}$, then stop and let ${\alpha _{inv}}:={\boldsymbol{\alpha }_{k}}$. Otherwise, update and go to the first step for the iteration again.
Algorithm A summarizes the proposed parameter estimation steps through the inverse problem approach which includes the details of the L-M method. In particular, the Armijo rule [2] in the third step, known as one of the inexact line search techniques, is imposed to ensure the objective function $\mathcal{G}$ has sufficient descent. Other rules and the related convergence theory can be found in [24].
3 An illustrative example
In this section, we use a heavy-tailed groundwater tracer concentration data to illustrate the methodology described in Section 2. The data comes from natural-gradient tracer tests conducted at the MacroDispersion Experimental (MADE) site at Columbus Air Force Base in northeastern Mississippi, precisely the MADE-2 tritium plume data [8]. The data consists of the maximum concentration measured in vertical slices perpendicular to the direction of plume travel, at day 27, day 132, day 224, and day 328 days after injection (see Figure 4 in [5]).
In [11], an fADE with constant drift and diffusion coefficient was fitted to the same data. With constant coefficient functions, the resulting solution process ${Y_{t}}$ in the SDE (1.2) itself is a stable process. The parameters of the fitted stable process were estimated using the least-square method.
Here we are considering an fADE that may have variable coefficients. For the fits presented in this section, the diffusion coefficient is still assumed to be location invariant, but we use a piecewise linear drift function of the form (2.24). The choice of ${x_{l}},\hspace{2.5pt}{x_{m}}$ and ${x_{r}}$ are subjective but should comply with the boundary conditions. ${x_{l}}=0$ (starting location), ${x_{r}}=300$ (beyond the observed location range) and ${x_{m}}=9.375$ (a rough mid-location of the observed range) were set for the current simulation. We used the fitted parameter values from the constant coefficient model [11] as initial values for the simulation; these parameter values are included in Appendix A. Then we used the data to estimate the linear drift and to obtain $p(x,t)$ that solves (2.1) by the methodology presented in Section 2. $c(x,t)=Kp(x,t)$ gives the fitted concentration that solves fADE associated with (2.1).
The concentrations fitted by the proposed method for the MADE site data at day 224 and day 328 are included here.
Fitted parameters for the MADE site data at day 224: ${a_{0}}=0.110$, ${a_{1}}=0.00032$, ${a_{2}}=0.0003$, ${a_{3}}=0.00019$, $b=0.1859783$, $\lambda =0.80$, $\gamma =0.9999$, $K=56778.24$.
Table 1.
Observed and fitted tracer concentrations for MADE data at day 224
| Day 224 | ||
| x = location, C = tracer concentration | ||
| x | Observed C | Fitted C |
| 2.1000 | 3378.0000 | 2871.7248 |
| 2.9000 | 1457.0000 | 3187.0093 |
| 3.6000 | 6494.0000 | 3529.6230 |
| 6.0000 | 1335.0000 | 4027.3618 |
| 6.6000 | 7705.0000 | 3321.0305 |
| 6.8000 | 2206.0000 | 3086.6518 |
| 6.9000 | 1291.0000 | 2969.4625 |
| 7.2000 | 4515.0000 | 2734.3773 |
| 8.0000 | 3598.0000 | 2361.1656 |
| 8.7000 | 2447.0000 | 2083.3103 |
| 9.7000 | 2831.0000 | 1809.0524 |
| 10.8000 | 2208.0000 | 1586.8821 |
| 12.7000 | 849.0000 | 1273.2312 |
| 13.5000 | 2213.0000 | 1169.6987 |
| 16.0000 | 1485.0000 | 935.2569 |
| 27.5000 | 443.0000 | 442.0004 |
| 31.8000 | 165.0000 | 356.0560 |
| 40.3000 | 291.0000 | 245.6054 |
| 48.8000 | 237.0000 | 178.6647 |
| 57.5000 | 76.0000 | 134.2397 |
| 66.3000 | 54.0000 | 103.7457 |
| 74.8000 | 137.0000 | 82.8644 |
| 83.4000 | 37.0000 | 67.3318 |
| 100.7000 | 28.0000 | 46.5417 |
| 117.0000 | 51.0000 | 34.4212 |
| 166.6000 | 18.0000 | 16.5498 |
| 183.5000 | 28.0000 | 13.9355 |
| 218.5000 | 9.0000 | 11.4791 |
| 268.7000 | 6.0000 | 10.1389 |
Fitted parameters for day 328 MADE site data: ${a_{0}}=0.105$, ${a_{1}}=0.00030$, ${a_{2}}=0.0005$, ${a_{3}}=0.00018$, $b=0.2233695$, $\lambda =0.79$, $\gamma =0.9999$, $K=37195.05$.
Table 2.
Observed and fitted tracer concentrations for MADE data at day 328
| Day 328 | ||
| x = location, C = tracer concentration | ||
| x | Observed C | Fitted C |
| 2.1000 | 1762.0000 | 2001.9586 |
| 3.6000 | 4157.0000 | 3264.5572 |
| 5.5000 | 2018.0000 | 1982.1889 |
| 6.0000 | 1295.0000 | 1653.7130 |
| 6.6000 | 2742.0000 | 1436.3243 |
| 6.8000 | 798.0000 | 1405.9525 |
| 6.9000 | 661.0000 | 1390.7665 |
| 7.2000 | 4721.0000 | 1323.8183 |
| 8.0000 | 2877.0000 | 1132.2689 |
| 8.7000 | 946.0000 | 1050.1787 |
| 9.7000 | 3140.0000 | 920.6741 |
| 10.8000 | 2075.0000 | 819.8101 |
| 12.7000 | 1636.0000 | 689.5556 |
| 13.5000 | 1569.0000 | 646.4434 |
| 16.0000 | 869.0000 | 536.8714 |
| 16.3000 | 1825.0000 | 525.8232 |
| 17.0000 | 286.0000 | 501.7569 |
| 18.0000 | 597.0000 | 470.2962 |
| 27.5000 | 127.0000 | 281.0554 |
| 31.8000 | 49.0000 | 231.1571 |
| 40.3000 | 108.0000 | 163.9576 |
| 48.8000 | 75.0000 | 121.3509 |
| 57.5000 | 95.0000 | 92.2235 |
| 66.3000 | 38.0000 | 71.8197 |
| 74.8000 | 126.0000 | 57.6513 |
| 83.4000 | 54.0000 | 47.0086 |
| 100.7000 | 19.0000 | 32.6278 |
| 117.0000 | 46.0000 | 24.1752 |
| 166.6000 | 25.0000 | 11.6326 |
| 183.5000 | 22.0000 | 9.4787 |
| 218.5000 | 11.0000 | 7.0942 |
| 268.7000 | 8.0000 | 6.0752 |
Comments:
-
1. The fit above uses piecewise linear drift velocity and location invariant diffusion coefficient as opposed to constant drift but same diffusion coefficient used in [11]. The fitted curves in Figures 2 and 4 show a better tail fit here compared to similar plots in [11]. The fractional order α for the fADE with constant coefficient was fitted as 1.09 and 1.05 for day 224 and day 328 respectively in [11]. Here the fitted fractional order $\alpha (=2-\lambda )$ are 1.2 and 1.19 for day 224 and day 328 respectively. Since the new fit improves on the fit in [11], this indicates that a better drift velocity fit may lead to an adjustment of the fractional order to improve the tail fit in the fADE. Ideally, the drift velocity function might be mechanically estimated from a properly designed field experiment while the diffusion coefficient and the fractional order can be estimated as described in this paper for an even better fADE fit.
-
2. Using simulation step 2, iid sample of size 1000 for ${Y_{t}}$ was generated at $t=\text{day}\hspace{2.5pt}224$ and $t=\text{day}\hspace{2.5pt}328$. Histograms for the generated samples in Figures 5 and 6 are fairly consistent with the observed concentration (or observed density) data. These figures indicate that although we had adopted a numerical approach to generate ${Y_{t}}$’s, the simulation method is efficient enough to follow the underlying probability distribution.
Fig. 5.
Histogram of 1000 ${Y_{t}}$ values generated by the proposed method at $t=224$. The red line shows the fitted pdf from (2.1) and blue dots show observed density = observed concentration /K (data estimates $K=56778.24$)
Fig. 6.
Histogram of 1000 ${Y_{t}}$ values generated by the proposed method at $t=328$. The red line shows the fitted pdf from (2.1) and blue dots show observed density = observed concentration /K (data estimates $K=37195.05$)
4 Discussion
The trapezoidal rule provides a very simple and quick numerical approximation of the CDF of ${Y_{t}}$ in (1.2) and for its inverse in the proposed simulation step 2. A more sophisticated numerical method like Simpson’s rule or other higher order Newton–Cotes formula or quadrature rules [23] can be used to obtain a more accurate ${\hat{F}_{{Y_{t}}}}$, but the procedure will be much more computationally taxing.
The method described here can be used to simulate a large number of observations (approximated) from ${Y_{t}}$. With the usual particle tracking insights, these simulated values can also be used to build empirical confidence intervals for the fitted concentration values.
Let us consider $P[{Y_{t}}\in (a,b)]={F_{{Y_{t}}}}(b)-{F_{{Y_{t}}}}(a)$. This is the probability that a randomly chosen tracer particle will be in a given interval $(a,b)$ at a time point t. For example, in groundwater pollution modeling this expression can be an important one that estimates the chance that the pollution will reach a certain area after a certain time. Using i.i.d. random simulations ${Y_{t}^{(1)}},{Y_{t}^{(2)}},\dots ,{Y_{t}^{(N)}}$ of ${Y_{t}}$, an empirical estimate for this probability can be given by $\hat{P}=\frac{1}{N}{\sum _{i=1}^{N}}I[{Y_{t}^{(i)}}\in (a,b)]$. Then by the central limit theorem $\hat{P}$ asymptotically follows the normal distribution with mean $P[{Y_{t}}\in (a,b)]$ and variance $\frac{P[{Y_{t}}\in (a,b)](1-P[{Y_{t}}\in (a,b)])}{N}$. Hence for large N’s, an asymptotic $(1-\alpha )100\% $ confidence interval for $P[{Y_{t}}\in (a,b)]$ can be given by
Now noting that density $p(y,t)$ can be approximated by $\frac{1}{b-a}P[{Y_{t}}\in (a,b)]$ with a small interval $(a,b)$, the simulated ${Y_{t}}$ values can be used to calculate the confidence intervals for the fitted densities. See [11] for detailed asymptotic confidence interval construction steps and related results associated with the empirical density function that is calculated using the generated ${Y_{t}}$’s.
Assuming that the observed concentrations are scaled versions of density $p(y,t)$ of ${Y_{t}}$, one can repeat simulation step 2 described in Section 2.1 without solving the forward equation to generate ${Y_{t}}$’s. However, the grid used for the simulation will be limited only to the fixed observed data location points. Since ${Y_{t}}$’s are simulated through numerical approximation the performance of this simulation depends on the grid spacing used. Thus the fitted fADE in simulation step 1 not only gives us the better understanding of the underlying process but also is essential for generating a good approximate process that resembles ${Y_{t}}$ by enabling the use of a finer grid for the numerical approximation. Further, modeling the underlying fADE process facilitates better prediction than mere extrapolation from an observed sample set.
On the other hand, an fADE can be fitted to the observed concentration by simulation step 1 only, without describing the underlying stochastic process. But the stochastic diffusion description is useful for understanding the mesoscopic flow and particle tracking methods. Further, generating the solution to the SDE associated with the fADE can be used to construct confidence intervals for the fADE fits that account for the error of estimation by a sample set.
In conclusion, the fADE and the associated α-stable SDE are essential tools to model heavy-tailed diffusion. However, the model fitting part may get complicated if we consider any scenario other than constant drift and diffusion coefficients. A numerical scheme for solving the fADE and the related SDE with linear drift is presented here. The general ideas used in the simulation steps can be replicated for a more complicated form of the drift. We plan to explore such models in our future work along with the convergence and the efficiency of the numerical approximations applied here.
