# Finite-difference form of the reaction-term in the solute transport equation

The partial differential equation is a combination of the diffusion plus convective trans­port equations and an adsorption sink. The equation for one-dimensional solute transport model is:

$$\frac{\partial C}{\partial t} = D\frac{\partial^ 2C}{\partial^ 2x}-v\frac{\partial C}{\partial x} - \frac{\rho}{\theta} \frac{\partial S}{\partial t}$$

where, C = solute concentration, D = dispersion coefficient, v = average pore-water velocity, x = distance from the inflow position, and t = time. Assuming the adsorption process is a first order reversible reaction, the rate of mass transfer to the adsorbed phase, $\frac{\partial S}{\partial t} = \frac{k_{A}\theta C}{\rho}-k_{D}S$; where, $k_{A}$ and $k_{D}$ are the adsorption (forward) and desorption (backward) rate coefficients (unit: 1/time), $\theta$ is the soil-water content by volume, and $\rho$ is the bulk density of the soil system.

The fully explicit finite-difference approximation for all except for the first order reversible reaction term can be written simply as (also, tested to work fine against exact solution):

$$C_{x,t} = C_{x,t-\Delta t} + \frac{D \Delta t}{\Delta x^2}(C|_{x+\Delta x, t-\Delta t}- 2C|_{x,t-\Delta t} + C|_{x-\Delta x,t-\Delta t}) - \frac{v\Delta t}{2 \Delta x}(C|_{x+\Delta x,t-\Delta t} - C|_{x-\Delta x, t-\Delta t})$$

I cannot seem to figure out how the above finite-difference approximation could be modified to incorporate the reaction-term defined above. Hint: Page#96-99 of the this book does provide a solution but I just cannot get my head around it. I'm supplying the best known articles for the analytical solution and numerical solution that I could find. Any help with reproducible example codes would be hihgly appreciated.

• If you supply BC, IC, and values for the parameters, I could try to obtain a solution. – Bill Greene Nov 10 '16 at 14:36
• What about $k_D$ and $k_A$? Also, $S$ really isn't a parameter but rather a dependent variable that will vary as a function of $x$ because it is a function of $C$ and also $t$. So is $S=.3$ an initial condition? – Bill Greene Nov 10 '16 at 15:31
• @BillGreene its the same conditions I indicated in my previous post: $x=0.075$ $m$, $D = 0.0033$ $m^2/day$, $v =0.66$ $m/day$, a constant initial concentration (left boundary) $C_{0} = 320$ $ppt$ (except that $C_{x} = 0$ when $t = 0$), and right boundary at an infinite/semi-infinite distance is 0. The $\Theta$, $\rho$, and $S$ parameters are: 0.3 cc/cc, 2.65 g/cc, and 1 ug/g, respectively. The constants $K_{A}$ and $K_{D}$ (unit: 1/day) could to be chosen as such they reproduce an analytical solution with retardation factor, R = 20, for example. – ToNoY Nov 10 '16 at 15:32
• no, S = 1 ug/g as initial condition, the porosity is 0.3 and bulk density of solid is 2.65 g/cc; sorry for mixing up units in cgs and metric system. – ToNoY Nov 10 '16 at 15:49
• @BillGreene shared two articles, see my updated question with two links – ToNoY Nov 13 '16 at 0:36

Typically you rewrite the first equation for $C$ to $$\frac{\partial C}{\partial t} = D\frac{\partial^ 2C}{\partial^ 2x}-v\frac{\partial C}{\partial x} - k_A C + \frac{\rho}{\theta} k_D S$$ using the second equation for $S$ that I repeat here $$\frac{\partial S}{\partial t} = \frac{k_{A}\theta}{\rho} C - k_{D} S \,.$$ If you consider the simplest (first order accurate in time) explicit numerical scheme for these two equations then using your notation it takes the form $$C_{x,t} = C_{x,t-\Delta t} + \frac{D \Delta t}{\Delta x^2}(C|_{x+\Delta x, t-\Delta t}- 2C|_{x,t-\Delta t} + C|_{x-\Delta x,t-\Delta t}) - \frac{v\Delta t}{2 \Delta x}(C|_{x+\Delta x,t-\Delta t} - C|_{x-\Delta x, t-\Delta t}) - \Delta t k_A C_{x,t-\Delta t} + \Delta t \frac{\rho}{\theta} k_D S_{x,t-\Delta t}$$ and $$S_{x,t}=S_{x,t-\Delta t} + \Delta t \frac{k_{A}\theta}{\rho} C_{x,t-\Delta t} - \Delta t k_{D} S_{x,t-\Delta t} \,.$$

P.S. Analogous equations are studied here:

De Smedt, F., Brevis, W., & Debels, P. (2005). Analytical solution for solute transport resulting from instantaneous injection in streams with transient storage. Journal of Hydrology, 315, 25–39.

http://users.clas.ufl.edu/jbmartin/website/Classes/Surface_Groundwater/Class%202/Smedt%20et%20al%20J%20Hydro%202005.pdf

Note that analytical solution is quite complicated, it involves so called special Bessel function and infinite sums. It should be analogous to the paper you quoted having analytical solution. If you plan only to compare numerical solution to available analytical solution, you might prefer only to compare published plots of analytical solution like Figure 1 in the above paper.

• can you explain a little why is it $C_{x,t-\Delta t}$ & $S_{x,t-\Delta t}$? Not involved with concentrations from any other nodes; e.g. $C_{x+\Delta x,t-\Delta t}$? The link I provided indicates a difference between two nodal concentrations divided by 2. Also, the $\theta$ term was with $C$, how come its with $S$ in your final equation? – ToNoY Nov 9 '16 at 3:21
• The link to the book you provided does not work for me, the pages are not visible. My first partial differential equation (PDE) for $C$ is obtained from your PDE simply by replacing $\frac{\partial S}{\partial t}$ there with the right hand side of ordinary differential equation (ODE) for $S$ that you have to multiply firstly by $\frac{\rho}{\theta}$. The second numerical scheme for $S_{x,t}$ is obtained by so called backward (explicit) Euler (time discretization) method that is the simplest possible one, you can use more sophisticated ones, but in fact you use it for the first PDE. – Peter Frolkovič Nov 9 '16 at 9:01
• Can you specify, please, which part in my answer is still unclear to you? Is the first differential equation for $C$ unclear? Is the second differential equation for $S$ unclear? Is the first numerical scheme for $C$ unclear? Is the second numerical scheme for S unclear? – Peter Frolkovič Nov 13 '16 at 18:27
• I'm not having any trouble understanding your solution but I also requested for a reproducible example of the scheme if possible. Once I can apply the whole scheme to a problem I'll accept the answer, but I'm working on it. – ToNoY Nov 13 '16 at 18:51
• OK, thanks for the answer. I have added a reference where also analytical solution to these equations is given. I think it is too complicated to implement only for the purpose of comparing numerical solution with the analytical one. You might prefer to reproduce plots of analytical solution that have been published for some particular parameters. – Peter Frolkovič Nov 13 '16 at 20:26