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é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.
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 ∂c∂t=−∂∂x[v(x)c]+∂2∂x2[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]:
∂c∂t=−∂∂x[v(x)c]+∂α−1∂xα−1[D(x)∂c∂x].
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 α∈(1,2) controls the tail of the flow. Here the fractional derivative of order α for any function f is defined as in [3] by:
dαf(x)dxα=−1Γ(1−α)∫0∞[f(x−u)−f(x)+uf′(x)]αy−1−αdu.
The negative fractional derivative is given by dαf(x)d(−x)α=dαg(x)dxα, 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 Xt1
Following the parametrization in [22] Xt∼Stable(α,β,μ,σ) where α is the index of stability, β is the skewness parameter, μ=0 is the location parameter and σ=1 is the scale parameter.
of the form:
dYt=a(Yt)dt+b(Yt)dXt.
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 Yt that solves SDE (1.2). It can be shown that Yt is a Markov process [1]. Let us denote the transition probability density function of Yt by py0(y,t), i.e. P(Yt∈A|Y0=y0)=∫Apy0(y,t)dy. For the initial distribution μ(u) the pdf of Yt is given by
p(y,t)=∫pu(y,t)μ(u)du,
e.g., for a ground water tracer concentration modeling it is reasonable to assume that all particles start at location y0 at time t=0 hence μ(u)=I(u=y0). If 1<α<2, it can be shown [10] that p(y,t) solves the forward equation:
∂p(y,t)∂t=−∂∂y[a(y)p(y,t)]+(1+β)2[−cos(πα2)]−1·∂α∂yα[bα(y)p(y,t)]+(1−β)2[−cos(πα2)]−1·∂α∂(−y)α[bα(y)p(y,t)].
Typically the coefficient functions a and b depend on t via Yt 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 Yt. However, without loss of generality, the functions can also be considered as a mapping from T×R→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:
∂C∂t=−∂∂x[v(x,t)C]+D(t)∂α∂xα[∂C∂x],
which essentially has the same form as the forward equation given in (1.4) with β=1, a(x,t)=v(x,t), [−cos(πα2)]−1bα(x,t)≡D(t). This shows that in a heavy-tailed plume that follows fADE (1.5), the position Yt of a randomly chosen tracer particle at time t solves SDE (1.2). Choosing β≠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)=Ktc(y,t), where Kt is a suitable scale parameter that adjusts the total mass for c(x,t) so that ∫Rp(y,t)dy=1 for each t.
In case v(·) and D(·) are constants, (1.2) reduces to Yt=at+bXt with a=v, b=[−cos(πα2)D]1/α and consequently, Yt 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 Yt 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 Yt using these fitted coefficient functions. The detailed methodology is described in Section 2. Section 3 provides an illustrative example. Section 4 includes concluding remarks.
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 Yt at a given time point t is given by p(y,t) and the cumulative distribution function (CDF) is defined by FYt(x)=∫−∞xp(y,t)dy. Since F(Yt)∼Uniform(0,1) [21], a random observation of Yt, say y, can be generated by FYt−1(u)=y where u is a random observation from 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 FYt and FYt−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(·) and D(·) 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 Xt 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 Yt using these p(t,y)’s.
Assuming the probability density p(y,t)=Ktc(y,t), we estimate Kt, α 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: GeneratingYt:
To generate Yt for a fixed t using an inverse CDF transformation, we present a straightforward scheme with a trapezoidal rule to approximate FYt and FYt−1. Approximate Yt’s can be simulated using p(y,t)’s from step 1 as follows:
Choose suitable grid y0<y1<⋯<yn and solve for p(yi,t) as mentioned in step 1. The endpoints y0 and yn can be chosen so that p(y,t) is small or negligible for y<y0 or y>yn{y_{n}}$]]>. The fitted coefficient and parameters from the observed data in step 1 can be used for any yi’s.
Approximate FYt(yi) by FˆYt(yi) for i=1,2,…,n, using the trapezoidal rule:
FˆYt(y1)=(y1−y0)p(y1,t)+p(y0,t)2,FˆYt(y2)=FˆYt(y1)+(y2−y1)p(y2,t)+p(y1,t)2,⋮FˆYt(yr)=FˆYt(yr−1)+(yr−yr−1)p(yr,t)+p(yr−1,t)2,⋮FˆYt(yn)=FˆYt(yn−1)+(yn−yn−1)p(yn,t)+p(yn−1,t)2.
For a randomly generated Uniform(0,1) observation u, find its placement in the grid, i.e. yr, yr+1 such that FˆYt(yr)≤u<FˆYt(yr+1), for some r∈{0,1,…,n}. A way of generating associated random observation of Yt or FYt−1(u) (approximated) as yˆ is described below:
yˆ=yr,ifu=yr;yr+u−FˆYt(yr)p(t,yr),ifyr<u<yr+1.
In case u<FˆYt(y0) set yˆ=y0 and if u≥FˆYt(yn) set y=yn.
To generate N i.i.d approximated random observations yˆ1,yˆ2,…,yˆN of Yt, start with randomly generated independent Uniform(0,1) observations u1,u2,…uN 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 FˆYt here and FˆYt→FYt as max0≤i≤n|yi+1−yi|→0 following the convergence property of the standard trapezoidal rule.
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−λ with 0<λ<1 (consistent with the construction described in [25]) as follows:
∂p∂t+∂∂xa(x)p−∂∂x[b(t)(γxlDx−λ+(1−γ)xDxr−λ)Dp]=0,xl<x<xr,p(xl,t)=0,p(xr,t)=0,0<t<T.
Essentially, compared to (1.4) the fractional derivative order α=2−λ, while γ=1+β2 indicates the relative weight of forward versus backward transition. xlDx−λ and xDxr−λ represent the left and right fractional integral operators defined as
xlDx−λf(x):=1Γ(β)∫xlx(x−s)λ−1f(s)ds,xDxr−λf(x):=1Γ(β)∫xxr(s−x)λ−1f(s)ds.
To comply with the boundary conditions, x=xl and x=xr are set to be the inflow and outflow boundaries respectively [25], while a(xl) and a(xr) are assumed to be non-negative.
The discretized finite volume scheme and the accumulation terms
Let us define a partition Pt on the time interval [0,T], and a partition Px on the space interval [xl,xr]Pt:0=t0<t1<⋯<tN=T,Px:xl=x0<x1<⋯<xI=xr,
with Δtn=tn−tn−1 for n=1,…,N and Δxi=xi−xi−1 for i=1,…,I. Let Px¯ be the dual partition of Px defined by xi−1/2=12(xi−1+xi) with Δxi+1/2=xi+1/2−xi−1/2 for i=1,…,I−1 as well as Δx1/2=x1/2−x0 and ΔxI+1/2=xI−xI−1/2.
Let y=r(θ;x¯,t¯) be a continuous and piecewise-smooth curve that passes through the point x¯ at time t¯, such that r(θ;xi−1/2,tn) and r(θ;xi+1/2,tn) do not meet each other during the time period [tn−1,tn]. We define a space-time control volume Ωin by extending the cell [xi−1/2,xi+1/2] along the curve r(t;x,tn) from t=tn to t=tn−1Ωin={(x,t):r(t;xi−1/2,tn)<x<r(t;xi+1/2,tn),tn−1<t<tn}.
Assuming the prism Ωin does not intersect the boundary x=xl and x=xr of the domain during the time period [tn−1,tn], let x∗=r(tn−1;x,tn) be the foot at time tn−1 of the curve with head x at time tn.
Integrating equation (2.1) over the control volume Ωin we get
∫Ωin∂p∂tdxdt+∫tn−1tn(a(x)p−b(t)∂1−λ∂x1−λp)|r(t;xi−1/2,tn)r(t;xi+1/2,tn)dt=0.
Without loss of generality, the accumulation term in (2.4) can be evaluated by assuming
xi−1/2∗<xi−1/2<xi+1/2∗<xi+1/2,
and accumulation term can be re-written as follows:
∫Ωin∂p∂tdxdt=∫xi−1/2∗xi−1/2[∫tn−1t(x;xi−1/2,tn)∂p∂tdt]dx+∫xi−1/2xi+1/2∗[∫tn−1tn∂p∂tdt]dx+∫xi+1/2∗xi+1/2[∫t(x;xi+1/2,tn)tn∂p∂tdt]dx.
Here the notation t(x;xi−1/2,tn) represents the time instant that the curve r(t;,xi−1/2,tn)=x. The notation t(x;xi+1/2,tn) is defined similarly.
A simple calculation of the second term on the right-hand side of (2.7) yields
∫xi−1/2xi+1/2∗[∫tn−1tn∂p∂tdt]dx=∫xi−1/2xi+1/2∗p(x,tn)dx−∫xi−1/2xi+1/2∗p(x,tn−1)dx.
The first and third terms on the right-hand side of (2.7) are integrated as follows
∫xi−1/2∗xi−1/2[∫tn−1t(x;xi−1/2,tn)∂p∂tdt]dx=∫tn−1tnp(r(t;xi−1/2,tn),t)∂r(t;xi−1/2,tn)∂tdt−∫xi−1/2∗xi−1/2p(x,tn−1)dx,
and
∫xi+1/2∗xi+1/2[∫t(x;xi+1/2,tn)tn∂p∂tdx]dt=∫xi+1/2∗xi+1/2p(x,tn)dx−∫tn−1tnp(r(t;xi+1/2,tn),t)∂r(t;xi+1/2,tn)∂tdt.
Incorporating equations (2.8)–(2.18) into (2.7) we get
∫Ωin∂p∂tdxdt=∫xi−1/2xi+1/2p(x,tn)dx−∫xi−1/2∗xi+1/2∗p(x,tn−1)dx+∫tn−1tnp(r(t;xi−1/2,tn),t)∂r(t;xi−1/2,tn)∂tdt−∫tn−1tnp(r(t;xi+1/2,tn),t)∂r(t;xi+1/2,tn)∂tdt.
The derivation shows that (2.11) does not depend on the assumption (2.4).
Now we just set the space-time boundaries r(t;xi±1/2,tn) of the control volume Ωin to be the characteristic curves, which are defined by the initial value problem of the ordinary differential equation
drdt=a(r,t),r(t;x,tn)|t=tn=x.
That is, because the characteristics r(t;xi±1/2,tn) 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 Ωin as follows:
∫xi−1/2xi+1/2p(x,tn)dx+∫tn−1tnF(p,b)(r(t;xi+1/2,tn),t)dt−∫tn−1tnF(p,b)(r(t;xi−1/2,tn),t)dt=∫xi−1/2∗xi+1/2∗p(x,tn−1)dx,
which can be approximated as
∫xi−1/2xi+1/2p(x,tn)dx+Δt[F(p,b)(xi+1/2,tn)−F(p,b)(xi−1/2,tn)]=∫xi−1/2∗xi+1/2∗p(x,tn−1)dx,
with the Lagrangian interface fluxes
F(p,b)(r(t;xi±1/2,tn),t):=−b∂1−λ∂x1−λp(r(t;xi±1/2,tn),t).
Note that the interface fluxes F(p,b)(r(t;xi±1/2,tn),t) are defined across the space-time boundary r(t;xi−1/2,tn) and r(t;xi+1/2,tn) of the space-time control volume Ωin.
To discretize equation (2.1), we further let {ϕi}i=1I be a set of hat functions such that ϕi(xi)=1 and ϕi(xj)=0 for j≠i. Hence the finite volume approximation ph to the true solution p can be expressed as
ph(x)=∑j=1I−1pjϕj(x).
and the finite volume scheme on the interval [xi−1/2,xi+1/2] for i=1,…,I−1 has the form
∫xi−1/2xi+1/2p(x,tn)dx+Δt∑j=1I−1pjnzi,j=∫xi−1/2∗xi+1/2∗p(x,tn−1)dx,
where, [zi,j]i,j=1I−1 is the coefficient matrix of the fractional term, and is given by
zi,j=b[(γ0Dxi−1/2−λ+(1−γ)x1−1/2D1−λ)Dϕj(xi−1/2,tn)−(γ0Dx+1/2−λ+(1−γ)x1+1/2D1−λ)Dϕj(xi+1/2,tn)].
Assuming the trial functions p(x,tn) are chosen to be piecewise linear functions on [xl,xr] with respect to the fixed spatial partition Px, we evaluate the accumulation term at time step tn in the finite volume scheme (2.16) analytically as
∫xi−1/2xi+1/2p(x,tn)dx=18[Δxi(pi−1n+3pin)+Δxi+1(3pin+pi+1n)].
In addition, we can compute xi−1/2∗ and xi+1/2∗ at time step tn−1 via a backward approximate characteristic tracking. Since the trial function p(x,tn−1) is also piecewise linear with respect to the fixed spatial partition Px at time tn−1, we can evaluate the accumulation term at time step tn−1 analytically. Because the accumulation term at time step tn−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.
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 ϕj have local support, (2.2) reveals that xlDx−λϕj and xDxr−λϕj have global support. Therefore, the stiffness matrix Z:=[zi,j]i,j=1I−1 is a full matrix, which requires O(N2) of memory to store. Numerical schemes for space-fractional differential equations were traditionally solved by the Gaussian type direct solvers that require O(N3) 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:
zi,i=1Γ(λ+1)2λbγhλ−1+1Γ(λ+1)b(1−γ)hλ−1(212λ−32λ)+1Γ(λ+1)bγhλ−1(212λ−32λ)+1Γ(λ+1)2λb(1−γ)hλ−1=1Γ(λ+1)bhλ−1(2−λ+212λ−32λ),
for 1≤i≤I−1.
Whereas, the sub-triangular entries of the matrix Z are given by:
zi,i−1=−1Γ(λ+1)bγhλ−1(212λ−32λ)−1Γ(λ+1)2λb(1−γ)hλ−1−1Γ(λ+1)bγhλ−1(52λ−232λ+12λ)=1Γ(λ+1)bhλ−1(332λγ−312λγ−52λγ−2−λ(1−γ)),
for 2≤i≤I−1, and
zi,j=1Γ(λ+1)bγhλ−1(i−j+12λ−2i−j−12λ+i−j−32λ)−1Γ(λ+1)bγhλ−1(i−j+32λ−2i−j+12λ+i−j−12λ)
for 3≤i≤I−1 and 1≤j≤i−2.
The super-triangular entries of matrix Z can be also derived as
zi,i+1=−1Γ(λ+1)b(1−γ)hλ−1(12λ−232λ+52λ)−1Γ(λ+1)2λbγhλ−1−1Γ(λ+1)b(1−γ)hλ−1(212λ−32λ)=1Γ(λ+1)bhλ−1(332λ(1−γ)−312λ(1−γ)−52λ(1−γ)−2−λγ),
for 1≤i≤I−2, and
zi,j=1Γ(λ+1)b(1−γ)hλ−1(2j−i+12λ−j−i−12λ−j−i+32λ)−1Γ(λ+1)b(1−γ)hλ−1(2j−i−12λ−j−i−32λ−j−i+12λ)
for 1≤j≤I−3 and j+2≤i≤I−1.
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,
a(x)=a0−a1x,ifxl≤x≤xm,a2−a3x,ifxm<x≤xr,
where parameters a0, a1, a2 and a3 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 p0 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 a0, a1, a2 and a3 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 a0 and a1. The other two parameters a2 and a3 can be estimated similarly. The parameter identification of a0 and a1 can be formulated as follows: let α:={a0,a1}, then to find αinv that satisfies
αinv=argminαG(α):=12∑i=1M[p(xi,t;α)−pˆ(xi,t)]2,
where, pˆ(xi,t)=Kˆtc(xi,t) for observed concentration c at location xi, time t.
In case the data is available for time points t1,t2,…,tR we rewrite equation (2.25) as
αinv=argminαG(α)=12∑k=1R∑i=1Mkwk[p(xtk,i,tk;α)−pˆ(xtki,tk)]2,
here xtk,i is the ith observed location at time tk and wk 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 G(α):
αk+1=αk−G′(αk)G″(αk),
where k represents the kth iteration. It is easy to check that (2.27) is equivalent to solve
αk+1=αk−(JkTJk)−1JkTrk,
where
Jk=(J1;J2;…;JR)T,withJi=∂p(xti,1,ti;α)∂α,…,∂p(xti,Mi,ti;α)∂α,
and
rk=(r1;r2;…;rR)T,withri=(ri,1,…,ri,Mi),ri,j=p(xti,j,ti;α)−pˆ(xti,j,ti),
for j=1,…,Mi, and i=1,2,…,R. Note that in practice, we always use the finite difference
p(xti,j,ti;α+δ)−p(xti,j,ti;α)δ
with a small enough δ to approximate the derivatives in (2.29).
However, the Newton method may fail to work because of JkTJk may be nearly zero. Therefore, the search direction dk:=−JkTrk/JkTJk 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
αk+1=αk−(JkTJk+ϱkI2)−1JkTrk,
where ϱk is a positive penalty parameter, and I2 is a 2×2 identity matrix. The method coincides with the Newton algorithm when ϱk=0; and it gives a step closer to the gradient descent direction when ϱk is large.
Algorithm A: Parameter Identification Algorithm:
Given the observation data g and an initial guess α0={a0,a1}0, choose ρ∈(0,1), σ∈(0,1/2), ϱ0>00$]]>, and δ small enough. For k=0,1,…,M, (let M=M1+M2+M3),
Solve the equation (2.1) corresponding to αk and αk+δ respectively to obtain p(·,·,αk) and p(·,·,αk+δ).
Compute Jk and rk, and update the search direction
dk:=−(JkTJk+ϱkI2)−1JkTrk.
Determine the search step ρm by Armijo rule:
G(αk+ρmdk)≤G(αk)+σρmdkJkTrk
where m is the smallest non-negative integer.
If |ρmdk|≤Tolerance level, then stop and let αinv:=αk. Otherwise, update
αk+1:=αk+ρmdk,ϱk+1:=ϱk/2,
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 G has sufficient descent. Other rules and the related convergence theory can be found in [24].
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 Yt 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 xl,xm and xr are subjective but should comply with the boundary conditions. xl=0 (starting location), xr=300 (beyond the observed location range) and xm=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:a0=0.110, a1=0.00032, a2=0.0003, a3=0.00019, b=0.1859783, λ=0.80, γ=0.9999, K=56778.24.
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:a0=0.105, a1=0.00030, a2=0.0005, a3=0.00018, b=0.2233695, λ=0.79, γ=0.9999, K=37195.05.
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
The fitted and observed concentration at day 224 for MADE site data
The fitted and observed concentration with log scale at day 224 for MADE site data
The fitted and observed concentration at day 328 for MADE site data
Comments:
Fitted and observed concentration with log scale at day 328 for MADE site data
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 α(=2−λ) 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.
Using simulation step 2, iid sample of size 1000 for Yt was generated at t=day224 and t=day328. 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 Yt’s, the simulation method is efficient enough to follow the underlying probability distribution.
Histogram of 1000 Yt 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)
Histogram of 1000 Yt 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)
Discussion
The trapezoidal rule provides a very simple and quick numerical approximation of the CDF of Yt 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 FˆYt, 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 Yt. 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[Yt∈(a,b)]=FYt(b)−FYt(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 Yt(1),Yt(2),…,Yt(N) of Yt, an empirical estimate for this probability can be given by Pˆ=1N∑i=1NI[Yt(i)∈(a,b)]. Then by the central limit theorem Pˆ asymptotically follows the normal distribution with mean P[Yt∈(a,b)] and variance P[Yt∈(a,b)](1−P[Yt∈(a,b)])N. Hence for large N’s, an asymptotic (1−α)100% confidence interval for P[Yt∈(a,b)] can be given by
Pˆ∓zα/2Pˆ(1−Pˆ)N
Now noting that density p(y,t) can be approximated by 1b−aP[Yt∈(a,b)] with a small interval (a,b), the simulated Yt 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 Yt’s.
Assuming that the observed concentrations are scaled versions of density p(y,t) of Yt, one can repeat simulation step 2 described in Section 2.1 without solving the forward equation to generate Yt’s. However, the grid used for the simulation will be limited only to the fixed observed data location points. Since Yt’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 Yt 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.
Appendix
In [11] the tracer trajectory Yt 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 α, skewness parameter β, location parameter v (same as the assumed constant drift velocity in (1.5)) and scale parameter σ, where σα=Dt|cos(πα/2)| with constant diffusion coefficient function D in (1.5).
The fitted values in [11] are:
Day 224 fitted values:σ=5.137167, μ=43.915430; β=0.99; α=1.0915; v=0.196051m/dayD=0.1859783mα/dayK=56778.24.
Day 328 fitted values:σ=5.380654; μ=74.677975; β=0.99; α=1.050998; v=0.2276768m/dayD=0.2233695mα/dayK=37195.05
These were used as initial values for the iterations in Section 3.
ReferencesApplebaum, D.: Lévy Processes and Stochastic Calculus. Cambridge University Press (2004). MR2072890. https://doi.org/10.1017/CBO9780511755323Armijo, L.: Minimization of functions having Lipschitz continuous partial derivatives. Pac. J. Math.16, 1–3 (1966). MR0191071Baeumer, B., Meerschaert, M.M.: Stochastic solutions for fractional Cauchy problems. Fract. Calc. Appl. Anal.4, 481–500 (2001). MR1874479Bear, J.: Dynamics of Fluids in Porous Media. Dover, Mineola, New York (1972)Benson, D.A., Schumer, R., Meerschaert, M.M., Wheatcraft, S.W.: Fractional dispersion, Lévy motion, and the MADE tracer tests. Transp. Porous Media42(1–2), 211–240 (2001). MR1948593. https://doi.org/10.1023/A:1006733002131Bhattacharya, R.N., Gupta, V.K.: A theoretical explanation of solute dispersion in saturated porous media at the Darcy scale. Water Resour. Res.19(4), 938–944 (1983)Bhattacharya, R.N., Gupta, V.K., Sposito, G.: On the stochastic foundation of the theory of water flow through unsaturated soil. Water Resour. Res.12(3) (1976)Boggs, J.M., Beard, L.M., Long, S.E., McGee, M.P.: Database for the second macro dispersion experiment (MADE-2). EPRI report TR-102072 (1993)Breiten, T., Simoncini, V., Stoll, M.: Low-rank solvers for fractional differential equations. Electron. Trans. Numer. Anal. (ETNA)45, 107–132 (2016). MR3498143Chakraborty, P.: A stochastic differential equation model with jumps for fractional advection and dispersion. J. Stat. Phys.136(3), 527–551 (2009). MR2529683. https://doi.org/10.1007/s10955-009-9794-1Chakraborty, P., Meerschaert, M.M., Lim, C.Y.: Parameter estimation for fractional transport: A particle tracking approach. Water Resour. Res.45(10), 10415 (2009)Chavent, G.: Nonlinear Least Squares for Inverse Problems: Theoretical Foundations and Step-by-Step Guide for Applications. Springer, Netherlands (2009). MR2554448del-Castillo-Negrete, D., Carreras, B.A., Lynch, V.E.: Fractional diffusion in plasma turbulence. Phys. Plasmas11, 3854 (2004)Ervin, V.J., Roop, J.P.: Variational formulation for the stationary fractional advection dispersion equation. Numer. Methods Partial Differ. Equ.22, 558–576 (2005). MR2212226. https://doi.org/10.1002/num.20112Fetter, C.W.: Applied Hydrogeology. Prentice Hall, New York (2001)Fu, H., Wang, H., Wang, Z.: POD/DEIM reduced-order modeling of time-fractional partial differential equations with applications in parameter identification. J. Sci. Comput.74, 220–243 (2018). MR3742877. https://doi.org/10.1007/s10915-017-0433-8Jia, J., Wang, H.: A preconditioned fast finite volume scheme for a fractional differential equation discretized on a locally refined composite mesh. J. Comput. Phys.299, 842–862 (2015). MR3384754. https://doi.org/10.1016/j.jcp.2015.06.028Liu, F., Anh, V., Turner, I.: Numerical solution of the space fractional Fokker-Plank equation. J. Comput. Appl. Math.166, 209–219 (2004). MR2057973. https://doi.org/10.1016/j.cam.2003.09.028Meerschaert, M., Tadjeran, C.: Finite difference approximations for fractional advection-dispersion flow equations. J. Comput. Appl. Math.172, 65–77 (2004). MR2091131. https://doi.org/10.1016/j.cam.2004.01.033Nocedal, J., Wright, S.J.: Numerical Optimization. Springer, New York, USA (2006). MR2244940Resnick, S.: A Probability Path. Birkhäuser, Boston (2003). MR1664717Samorodnitsky, G., Taqqu, M.S.: Stable Non-Gaussian Random Processes, Stochastic Models with Infinite Variance. Chapman & Hall, New York (1994). MR1280932Scarborough, J.B.: Numerical Mathematical Analysis. Johns Hopkins Press (1966). MR0198651Sun, W., Yuan, Y.: Optimization Theory and Methods: Nonlinear Programming. Springer, New York, USA (2006). MR2232297Wang, H., Al-Lawatia, M.: A locally conservative Eulerian-Lagrangian control-volume method for transient advection-diffusion equations. Numer. Methods Partial Differ. Equ.22, 577–599 (2006). MR2212227. https://doi.org/10.1002/num.20106Wang, H., Du, N.: A superfast-preconditioned iterative method for steady-state space-fractional diffusion equations. J. Comput. Phys.240, 49–57 (2013). MR3039244. https://doi.org/10.1016/j.jcp.2012.07.045Zhang, Y., Benson, D., Meerschaert, M.M., LaBolle, E.: Space-fractional advection-dispersion equations with variable parameters: Diverse formulas, numerical solutions, and application to the MADE-site data. Water Resour. Res.43(5), 05439 (2007)