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

Question

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

Answer 1

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.