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