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

Question

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.

Answer 1

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.