This material is based upon work supported by the National Science Foundation under Grant Nos. 9980489 and 0139943. Latest update: 6/3/03

Sure, I'll code that right up for you ...

1. BACKGROUND

where P=P(x,t) is probability (i.e., concentration), B is the "capacity coefficient" decsribing the ratio of immobile to mobile porosity, g (between 0 and 1) is the fractional order of differentiation that describes the heavy-tailed waiting times of particles in that immobile phase, v is the mean velocity, D_F is the Fickian dispersion coefficient, D is the scalar anomalous dispersion coefficient, and the term that it is attached to is a matrix-order Laplacian. The order of differentiation in the fractional Laplacian varies with direction and is given by the inverse of the "scaling matrix" H, whose the eigenvalues are similar to the Hurst coefficients. The eigenvectors of H are the directions of priciple plume growth (not necessarily orthogonal). The mixing measure M(dq) describes the propensity for particle motions in any direction q on the unit sphere in d-dimensions. For example, in 1-D, particles can only go forward or backward, so dq = ± 1, and M(+1) = 0.5(1+beta), and M(-1) = 0.5(1-beta), and beta between -1 and 1 is the skewness in 1-D. For a Brownian motion, the fractional Laplacian term is superfluous, but it is intructive to see that it recovers classical diffusion. In this case, the scaling (Hurst) matrix H is diag(1/2), so H^-1 = 2I, recovering the Laplacian operator. An explanation of the spatial fractional operators can be found in:

R. Schumer, D.A. Benson, M.M. Meerschaert, and B. Baeumer, Multiscaling fractional advection-dispersion equations and their solutions (460 Kb pdf), Water Resources Research 39(1), 1022, doi:10.1029/2001WR001229, 2003.

An explanation of the time operators and the link to a fractal Mobile/Immobile model is given by Rina Schumer in her Ph.D. thesis.

2. STOCHASTIC UNDERPINNINGS

It suffices to solve the problem with only the first time derivative part, and then subordinate that solution accoding to the fractional time derivative part. Consequently, we have no all-encompassing solution to the last equation. Instead we solve it in two parts. A note to mathematicians: The following programs generate the probability density functions of 1 and 2-D stables and operator stables via Fourier transform. The tails of the operator stable are not assurate below about 4 decimal places.

3. SIMPLIFICATIONS (with software)

3.1 The spatial process (Cauchy Problems)

1-D, with automated fitting of parameters

In hydrologic applications, we suggest setting beta = 1, which models forward particle motions only. With this fortran program, you give intial guesses of mean velocity (v), dispersion coefficient (D), the skewness (beta), the order of the fractional space derivative (a between 0 and including 2), and/or initial mass (m_0). You may hold any of these constant. The program fits your data, either a set of C(t) at given x, or C(x) at a given time. We provide the core of the MADE site tritium plume as an example. Note that we wrote this fitting program after eyeball fitting in the paper Fractional dispersion, Lévy motion, and the MADE tracer tests (222 kB .pdf) Transp. Por. Media 42(1/2), 211-240, 2001.

2-D, with any orientation of eigenvectors of the plume growth process H

This Mathcad sheet uses Fourier transform to invert the known characteristic function of an operator stable Law. In other words, we solve (1) with any M(dq) and matrix H, with B=0 (no fractional time). The zip file contains two Mathcad sheets, one for orthogonal eigenvectors (granular media) and one for non-orthogonal growth eigenvectors (fractured media). Instructions are contained at the bottom of the sheets. If you do not have Mathcad (a few $100), you can at least be happy that we didn't use Mathematica, which goes for ~$1,000 a pop.

3.2 Subordination (the time process)

Combined with 1-D transport

This Mathcad sheet solves a simple 1-D transport equation with heavy-tailed "sorption" to an immobile phase. It first solves eq. (2), with either Fickian (2nd-order) or fractional spatial term. Then it subordinates the solution. We provide a subset of the stream dye experiment from Haggerty, R., S.M. Wondzell, M.A. Johnson, Power-law residence time distribution in the hyporheic zone of a 2nd-order mountain stream, Geophys. Res. Lett. 29(10), 2002. We thank those authors for allowing us to have and share the data (go get the data!).   Thanks also to Boris Baeumer (U. Otago, New Zealand) for Mathcad coding! Latest update May 2003 Read about the theory and application in Fractal Mobile/Immobile Transport (Rina Schumer, David A. Benson, Mark M. Meerschaert, and Boris Baeumer) (1057 Kb pdf) Water Resources Research 39(10)1296, doi:10.1029/2003WR02141, 2003.

General Subordination

This Mathcad sheet takes your solution to a Cauchy problem a(t) at some point x, and transforms it according to the subordination integral (adding the effects of havy tailed waiting times). You only specify B and g in eq. (1). It assumes that your time series is zero for time later than your final data point.  Latest update 3 June 2003.