# Unrealistic solution for advection-diffusion-reaction PDE with heterogeneous porous media

About the code: I have a code which simulates concentration from advection-diffusion-reaction PDE in 2D space (X,Y) with time. The solution is obtained using fully implicit finite-difference method and includes the capability to simulate a media with spatially varying permeability and reaction constant (through upwinding by harmonic mean). I have been able to test the code for homogeneous media and it works fine. It's a solution for the following equation: \begin{align} \frac{\partial C}{\partial t} + \nabla.\left(v C - D\nabla{C} \right)= \alpha C \end{align}

Issue: The code gives realistic solution for a medium with uniform layers of permeability as shown below. However, once the medium starts to get more heterogeneous it starts to throw unrealistic results (i.e. negative concentration and large fluctuations with time). Any guesses why I am getting unrealistic results for heterogeneous media (with values in the similar range as uniform layered media)?

UPDATE: I think I have figured the issue for getting negative values and sharp fluctuations in concentration. I think the problem is not with the numerical model, but with the physical values. Even though I was already taking realistic parameter values for which I was getting unrealistic concentration, I had to further adjust them (decrease the convection velocity) to get realistic results

• could be related to which finite difference scheme is being used(forward,backward,difference).could also be due to not using weighted values – shrey Nov 22 '15 at 5:31
• Forward difference in time and space. May I know which weights do you mean? – user3153 Nov 22 '15 at 6:32
• Not an expert, but I wonder if the problem is in your solver for the finite difference. IIRC, the transport equations can be subject to a large amount of numerical dispersion. What happens if you use a finer grid? – haresfur Nov 22 '15 at 21:57
• haresfur: I have tried finer grids but that didn't help either, though it changed the results to some extent but the problem of sharp oscillations still existed even with finer grids. – user3153 Nov 22 '15 at 22:04

## 1 Answer

It's not clear exactly what is being modelled here, but it seems to me that there are two ways in which the concentration can 'go negative'. Firstly, the rate of change of concentration can be massive, in which case see what happens when modelling with much smaller time steps. Or, the diffusion term substantially exceeds the advection term, which is physically unlikely. However, you are modelling a large (more than 10-fold) anisotropy, in which almost all the flow will be in the x, horizontal direction - it is almost a one-dimensional problem in fact. So how can the advection be less than the diffusion? In the case of a constant concentration source, such as a saline water body, the negative results don't make sense, in which case the calculations have gone unstable, probably from using too large a time step. In the case of a pulsed source, as the concentration peak passes (along the x direction), the advected concentration rises and falls, such that the lateral (y axis) diffusion increases, then decreases, i.e. the concentration gradient changes direction, becoming negative.

As if this wasn't hard enough, in real life applications the 'sink' terms probably far exceed any back-diffusive effects.

• Gordon: I have implemented fully implicit scheme in which the accuracy of results are independent from the effect of time step size. Additionally, even if diffusion is zero (in most cases it's almost negligible compared to advective flux), I still have this problem. – user3153 Nov 22 '15 at 18:19
• Then I'm a bit puzzled also. Do you have the source code, or is this an off the shelf package? – Gordon Stanger Nov 22 '15 at 19:24
• I wrote it myself...so I have the source code, Gordon. – user3153 Nov 22 '15 at 20:55
• In what language? Fortran I can read. C, I can't – Gordon Stanger Nov 23 '15 at 7:04
• In python, Gordon – user3153 Nov 23 '15 at 13:12