Stable L\'{e}vy diffusion and related model fitting

TL;DR

This paper introduces a numerical framework for solving stable Lévy stochastic differential equations related to fractional advection-dispersion equations, enabling modeling of heavy-tailed flows with unknown coefficients from observed data.

Contribution

It proposes a novel numerical method combining finite volume schemes and Levenberg--Marquardt regularization to solve and fit stable SDEs for complex flow data.

Findings

Numerical method effectively solves stable SDEs in general scenarios.

Levenberg--Marquardt regularization estimates unknown advection and dispersion functions.

Framework accurately models heavy-tailed flow phenomena.

Abstract

A fractional advection-dispersion equation (fADE) has been advocated for heavy-tailed flows where the usual Brownian diffusion models fail. A stochastic differential equation (SDE) driven by a stable L\'{e}vy process gives a forward equation that matches the space-fractional advection-dispersion equation and thus gives the stochastic framework of particle tracking for heavy-tailed flows. For constant advection and dispersion coefficient functions, the solution to such SDE itself is a stable process and can be derived easily by least square parameter fitting from the observed flow concentration data. However, in a more generalized scenario, a closed form for the solution to a stable SDE may not exist. We propose a numerical method for solving/generating a stable SDE in a general set-up. The method incorporates a discretized finite volume scheme with the characteristic line to solve the fADE or the forward equation for the Markov process that solves the stable SDE. Then we use a numerical scheme to generate the solution to the governing SDE using the fADE solution. Also, often the functional form of the advection or dispersion coefficients are not known for a given plume concentration data to start with. We use a Levenberg--Marquardt (L-M) regularization method to estimate advection and dispersion coefficient function from the observed data (we present the case for a linear advection) and proceed with the SDE solution construction described above.

Tables

Table 1. 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

Table 2. 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

Equations

$$\dfrac{\partial c}{\partial t}=-\dfrac{\partial}{\partial x}\bigl{[}v(x)c\bigr{]}+\dfrac{\partial^{\alpha-1}}{\partial x^{\alpha-1}}\biggl{[}D(x)\dfrac{\partial c}{\partial x}\biggr{]}.$$

$$\frac{d^{\alpha}f(x)}{dx^{\alpha}}=\frac{-1}{\varGamma(1-\alpha)}\int\limits^{\infty}_{0}\bigl{[}f(x-u)-f(x)+uf^{\prime}(x)\bigr{]}\alpha y^{-1-\alpha}du.$$

$$d{Y}_t=a(Y_t)dt+b(Y_t)d{X}_t.$$

$$p(y,t)=\int p_u(y,t)\mu (u)du,$$

$$\displaystyle\frac{\partial p(y,t)}{\partial t}=-\frac{\partial}{\partial y}\bigl{[}a(y)p(y,t)\bigr{]}$$

$$\displaystyle+\frac{(1-\beta)}{2}\biggl{[}\!-\cos\biggl{(}\frac{\pi\alpha}{2}\biggr{)}\!\biggr{]}^{-1}\!\!\cdot\frac{\partial^{\alpha}}{\partial(-y)^{\alpha}}\bigl{[}b^{\alpha}(y)p(y,t)\bigr{]}.$$

$$\dfrac{\partial C}{\partial t}=-\dfrac{\partial}{\partial x}\bigl{[}v(x,t)C\bigr{]}+D(t)\dfrac{\partial^{\alpha}}{\partial x^{\alpha}}\biggl{[}\dfrac{\partial C}{\partial x}\biggr{]},$$

$${\hat{F}}_{Y_t}(y_1)$$

$${\hat{F}}_{Y_t}(y_2)$$

$$⋮$$

$${\hat{F}}_{Y_t}(y_r)$$

$$⋮$$

$${\hat{F}}_{Y_t}(y_n)$$

$$\hat{y}=\left\{\begin{array}[]{ll}y_{r},&\hbox{if }\ u=y_{r};\\ y_{r}+\dfrac{u-\hat{F}_{Y_{t}}(y_{r})}{p(t,y_{r})},&\hbox{if }\ y_{r}<u<y_{r+1}.\end{array}\right.$$

$$\displaystyle\frac{\partial p}{\partial t}+\frac{\partial}{\partial x}a(x)p-\frac{\partial}{\partial x}\bigl{[}b(t)\bigl{(}\gamma{}_{x_{l}}D_{x}^{-\lambda}+(1-\gamma){}_{x}D_{x_{r}}^{-\lambda}\bigr{)}Dp\bigr{]}=0,\quad x_{l}<x<x_{r},$$

$$p(x_l,t)=0,p(x_r,t)=0,0<t<T.$$

$${}_{x_l}{D}_x^{-\lambda }f(x):$$

$${}_x{D}_{x_r}^{-\lambda }f(x):$$

$${\mathcal{P}}_t:0$$

$${\mathcal{P}}_x:x_l$$

$${\varOmega}_{i}^{n}=\bigl{\{}(x,t):r(t;x_{i-1/2},t_{n})<x<r(t;x_{i+1/2},t_{n}),\quad t_{n-1}<t<t_{n}\bigr{\}}.$$

$$\int_{{\varOmega}_{i}^{n}}\frac{\partial p}{\partial t}dxdt+\int_{t_{n-1}}^{t_{n}}\Bigl{(}a(x)p-b(t)\frac{\partial^{1-\lambda}}{\partial x^{1-\lambda}}p\Bigr{)}\Big{|}_{r(t;x_{i-1/2},t_{n})}^{r(t;x_{i+1/2},t_{n})}dt=0.$$

$$x_{i-1/2}^{\ast }<x_{i-1/2}<x_{i+1/2}^{\ast }<x_{i+1/2},$$

$${\int }_{{\Omega }_i^n}\frac{\partial p}{\partial t}dxdt$$

$$\displaystyle\quad+\int_{x_{i-1/2}}^{x_{i+1/2}^{*}}\Biggl{[}\int_{t_{n-1}}^{t_{n}}\frac{\partial p}{\partial t}\>dt\Biggr{]}\>dx+\int_{x_{i+1/2}^{*}}^{x_{i+1/2}}\Biggl{[}\int_{t(x;x_{i+1/2},t_{n})}^{t_{n}}\frac{\partial p}{\partial t}\>dt\Biggr{]}\>dx.$$

$$\int_{x_{i-1/2}}^{x_{i+1/2}^{*}}\Biggl{[}\int_{t_{n-1}}^{t_{n}}\frac{\partial p}{\partial t}\>dt\Biggr{]}dx=\int_{x_{i-1/2}}^{x_{i+1/2}^{*}}p(x,t_{n})dx-\int_{x_{i-1/2}}^{x_{i+1/2}^{*}}p(x,t_{n-1})dx.$$

$$\displaystyle\int_{x_{i-1/2}^{*}}^{x_{i-1/2}}\Biggl{[}\int_{t_{n-1}}^{t(x;x_{i-1/2},t_{n})}\frac{\partial p}{\partial t}\>dt\Biggr{]}dx$$

$$\displaystyle\quad=\int_{t_{n-1}}^{t_{n}}p\bigl{(}r(t;x_{i-1/2},t_{n}),t\bigr{)}\frac{\partial r(t;x_{i-1/2},t_{n})}{\partial t}\>dt-\int_{x_{i-1/2}^{*}}^{x_{i-1/2}}p(x,t_{n-1})\>dx,$$

$$\displaystyle\int_{x_{i+1/2}^{*}}^{x_{i+1/2}}\Biggl{[}\int_{t(x;x_{i+1/2},t_{n})}^{t_{n}}\frac{\partial p}{\partial t}\>dx\Biggr{]}dt$$

$$\displaystyle\quad=\int_{x_{i+1/2}^{*}}^{x_{i+1/2}}p(x,t_{n})\>dx-\int_{t_{n-1}}^{t_{n}}p\bigl{(}r(t;x_{i+1/2},t_{n}),t\bigr{)}\frac{\partial r(t;x_{i+1/2},t_{n})}{\partial t}~{}dt.$$

$${\int }_{{\Omega }_i^n}\frac{\partial p}{\partial t}dxdt$$

$$\displaystyle\quad+\int_{t_{n-1}}^{t_{n}}p\bigl{(}r(t;x_{i-1/2},t_{n}),t\bigr{)}\frac{\partial r(t;x_{i-1/2},t_{n})}{\partial t}\>dt$$

$$\displaystyle\quad-\int_{t_{n-1}}^{t_{n}}p\bigl{(}r(t;x_{i+1/2},t_{n}),t\bigr{)}\frac{\partial r(t;x_{i+1/2},t_{n})}{\partial t}\>dt.$$

$$\frac{dr}{dt}=a(r,t),r(t;x,t_n){|}_{t=t_n}=x.$$

$$\displaystyle\int_{x_{i-1/2}}^{x_{i+1/2}}p(x,t_{n})\>dx+\int_{t_{n-1}}^{t_{n}}F(p,b)\bigl{(}r(t;x_{i+1/2},t_{n}),t\bigr{)}\>dt$$

$$\displaystyle\qquad-\int_{t_{n-1}}^{t_{n}}F(p,b)\bigl{(}r(t;x_{i-1/2},t_{n}),t\bigr{)}\>dt$$

$$={\int }_{x_{i-1/2}^{\ast }}^{x_{i+1/2}^{\ast }}p(x,t_{n-1})dx,$$

$$\displaystyle\int_{x_{i-1/2}}^{x_{i+1/2}}p(x,t_{n})\>dx+\Delta t\bigl{[}F(p,b)(x_{i+1/2},t_{n})-F(p,b)(x_{i-1/2},t_{n})\bigr{]}$$

Full text

Stable Lévy Diffusion,

fractional diffusion,

fractional advection-dispersion,

heavy-tailed particle tracking,

\aid

VMSTA106

\startlocaldefs

\urlstylerm \endlocaldefs

\pretitle

Research Article

[type=corresp,id=cor1]Corresponding author.

\publishedonline\sday

9 7 2018

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]:

[TABLE]

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 $\alpha$ for any function $f$ is defined as in [3] by:

[TABLE]

The negative fractional derivative is given by $\frac{d^{\alpha}f(x)}{d(-x)^{\alpha}}=\frac{d^{\alpha}g(x)}{dx^{\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 $\alpha$-stable Lévy process $X_{t}$ 111Following the parametrization in [22] $X_{t}\sim\mathit{Stable}(\alpha,\beta,\mu,\sigma)$ where $\alpha$ is the index of stability, $\beta$ is the skewness parameter, $\mu=0$ is the location parameter and $\sigma=1$ is the scale parameter. of the form:

[TABLE]

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

[TABLE]

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:

[TABLE]

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}\rightarrow\mathbb{R}$ as $a(y,t)$ and $b(y,t)$ and the associated forward equation will remain the same.

Thus for a location invariant $D$, the fADE in (1.1) can be written as:

[TABLE]

which essentially has the same form as the forward equation given in (1) 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\neq 1$ will allow the particle transition in forward or backward direction.

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+bX_{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^{x}_{-\infty}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 $\alpha$-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) 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}$, $\alpha$ and $\beta$ 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). In case the coefficient functions in (1) 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:

1. (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. 2. (b)

Approximate $F_{Y_{t}}(y_{i})$ by $\hat{F}_{Y_{t}}(y_{i})$ for $i=1,2,\ldots,n$, using the trapezoidal rule:

[TABLE] 3. (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})\leq u<\hat{F}_{Y_{t}}(y_{r+1})$, for some $r\in\{0,1,\ldots,n\}$. A way of generating associated random observation of $Y_{t}$ or $F^{-1}_{Y_{t}}(u)$ (approximated) as $\hat{y}$ is described below:

[TABLE]

In case $u<\hat{F}_{Y_{t}}(y_{0})$ set $\hat{y}=y_{0}$ and if $u\geq\hat{F}_{Y_{t}}(y_{n})$ set $y=y_{n}$. 4. (d)

To generate $N$ i.i.d approximated random observations $\hat{y}_{1},\hat{y}_{2},\ldots,\hat{y}_{N}$ of $Y_{t}$, start with randomly generated independent $\mathit{Uniform}(0,1)$ observations $u_{1},u_{2},\ldots 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}}\rightarrow F_{Y_{t}}$ as $\max\limits_{0\leq i\leq n}|y_{i+1}-y_{i}|\rightarrow 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). 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:

[TABLE]

Essentially, compared to (1) 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

[TABLE]

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.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}]$

[TABLE]

with $\Delta t_{n}=t_{n}-t_{n-1}$ for $n=1,\ldots,N$ and $\Delta x_{i}=x_{i}-x_{i-1}$ for $i=1,\ldots,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,\ldots,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}$.

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}$

[TABLE]

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^{*}=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}$.

Integrating equation (2.1) over the control volume ${\varOmega}_{i}^{n}$ we get

[TABLE]

Without loss of generality, the accumulation term in (2.4) can be evaluated by assuming

[TABLE]

and accumulation term can be re-written as follows:

[TABLE]

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.

A simple calculation of the second term on the right-hand side of (2.2.1) yields

[TABLE]

The first and third terms on the right-hand side of (2.2.1) are integrated as follows

[TABLE]

and

[TABLE]

Incorporating equations (2.8)–(2.18) into (2.2.1) we get

[TABLE]

The derivation shows that (2.11) does not depend on the assumption (2.4).

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

[TABLE]

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:

[TABLE]

which can be approximated as

[TABLE]

with the Lagrangian interface fluxes

[TABLE]

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\neq i$. Hence the finite volume approximation $p_{h}$ to the true solution $p$ can be expressed as

[TABLE]

and the finite volume scheme on the interval $[x_{i-1/2},x_{i+1/2}]$ for $i=1,\ldots,I-1$ has the form

[TABLE]

where, $[z_{i,j}]_{i,j=1}^{I-1}$ is the coefficient matrix of the fractional term, and is given by

[TABLE]

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

[TABLE]

In addition, we can compute $x_{i-1/2}^{*}$ and $x_{i+1/2}^{*}$ 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:

[TABLE]

for $1\leq i\leq I-1$.

Whereas, the sub-triangular entries of the matrix $Z$ are given by:

[TABLE]

for $2\leq i\leq I-1$, and

[TABLE]

for $3\leq i\leq I-1$ and $1\leq j\leq i-2$.

The super-triangular entries of matrix $Z$ can be also derived as

[TABLE]

for $1\leq i\leq I-2$, and

[TABLE]

for $1\leq j\leq I-3$ and $j+2\leq i\leq I-1$.

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$,

[TABLE]

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 $\mathbf{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).

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

[TABLE]

where, $\hat{p}(x_{i},t)=\hat{K}_{t}c(x_{i},t)$ for observed concentration $c$ at location $x_{i}$, time $t$.

In case the data is available for time points $t_{1},t_{2},\ldots,t_{R}$ we rewrite equation (2.25) as

[TABLE]

here $x_{t_{k},i}$ is the $i$th 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$).

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})$:

[TABLE]

where $k$ represents the $k$th iteration. It is easy to check that (2.27) is equivalent to solve

[TABLE]

where

[TABLE]

and

[TABLE]

for $j=1,\ldots,M_{i}$, and $i=1,2,\ldots,R$. Note that in practice, we always use the finite difference

[TABLE]

with a small enough $\delta$ 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

[TABLE]

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.

Algorithm A: Parameter Identification Algorithm**:**

Given the observation data $\mathbf{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 $\delta$ small enough. For $k=0,1,\ldots,M$, (let $M=M_{1}+M_{2}+M_{3}$),

1. (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)$. 2. (2)

Compute $\mathbf{J}_{k}$ and $\mathbf{r}_{k}$, and update the search direction

[TABLE] 3. (3)

Determine the search step $\rho^{m}$ by Armijo rule:

[TABLE]

where $m$ is the smallest non-negative integer. 4. (4)

If $|\rho^{m}\mathbf{d}_{k}|\leq\text{Tolerance level}$, then stop and let $\alpha_{inv}:=\boldsymbol{\alpha}_{k}$. Otherwise, update

[TABLE]

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},\ 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$.

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$.

Comments:

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 $\alpha$ 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. 2.

Using simulation step 2, iid sample of size 1000 for $Y_{t}$ was generated at $t=\mbox{day }224$ and $t=\mbox{day }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.

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^{(1)}_{t},Y^{(2)}_{t},\ldots,Y^{(N)}_{t}$ of $Y_{t}$, an empirical estimate for this probability can be given by $\hat{P}=\frac{1}{N}\sum^{N}_{i=1}I[Y^{(i)}_{t}\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 $\dfrac{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

[TABLE]

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 $\alpha$-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.

Appendix A Appendix

In [11] the tracer trajectory $Y_{t}$ of the MADE site data was assumed to follow an SDE of the form (1.2) with constant parameters and hence itself was a stable process with index of stability $\alpha$, skewness parameter $\beta$, location parameter $v$ (same as the assumed constant drift velocity in (1.5)) and scale parameter $\sigma$, where $\sigma^{\alpha}=Dt|\cos(\pi\alpha/2)|$ with constant diffusion coefficient function $D$ in (1.5).

The fitted values in [11] are:

Day 224 fitted values: $\sigma=5.137167$, $\mu=43.915430$; $\beta=0.99$; $\alpha=1.0915$; $v=0.196051m/\textit{day}$ $D=0.1859783m^{\alpha}/\textit{day}$ $K=56778.24$.

Day 328 fitted values: $\sigma=5.380654$; $\mu=74.677975$; $\beta=0.99$; $\alpha=1.050998$; $v=0.2276768m/\textit{day}$ $D=0.2233695m^{\alpha}/\textit{day}$ $K=37195.05$

These were used as initial values for the iterations in Section 3.