Negative concentrations may be considered as a „self-diagnostic“ indicating weaknesses in the model setup.
There are several ways on dealing with negative concentrations. These different ways may be categorized in adaption and smoothing techniques. Whether adaption or smoothing techniques are preferred depends very much on the study area, individual model settings and also on the question you want to answer.
During the first attempts, I suggest prioritizing adaption techniques rather than smoothing techniques. If all attempts fail, you could try to use smoothing techniques.
Here are some (general) questions I would try to answer before I switch to smoothing techniques:
1. Is the mesh fine enough? A fine mesh is required especially at sharp interfaces of primary variables (e.g. heads or freshwater-saltwater interface) and in areas where sharp contrasts in physical rock properties characterize the hydrogeological setting.
2. Does the mesh contain bad shaped elements? Try to use the maximum interior angle of triangles instead of the Delaunay criterion violations. You may improve the mesh using mesh smoothing techniques from the Mesh Geometry Toolbar.
3. What temporal discretization are you using? I guess you use an automatic predictor-corrector time stepping scheme. That’s fine. But sometimes additional constraints by means of a growth factor between subsequent time-steps or/and of a maximum time-step size may improve your model.
4. Fully implicit schemes are generally more stable than (semi-)implicit schemes. Accordingly, you may switch from a second order-accurate semi-implicit scheme (AB/TR) to a fully implicit first-order accurate (FE/BE) scheme.
5. What error norm are you using? The maximum error norm represents a more stringent criterion compared to integral error norms given the fact that the maximum error norm focuses more on local effects within the model. In this regard you may also change the error tolerance.
Finally, if all these settings still do not give rise to the desired result or if (for instance) the mesh would contain too many elements, which exceeds the handling of your computational resources, you may think about adopting Upwinding to stabilize the numerical solution artificially. Upwinding techniques smooth steep gradients of computational findings by adding (artificial) numerical dispersion. In other words, you change the physics of your system in order to reach a numerically stable solution. Therefore, I always would be cautious when dealing with Upwinding.
You may also think about accepting negative concentrations to a certain degree. Criteria to assess whether or not a concentration is acceptable could be the magnitude of the concentration or the spatial location where the negative concentration occurs.
