I am running a 2D vertical, transient model of flow and mass transport aiming to simulate the Henry problem under different conditions. I run the simulation until the system reaches pseudo steady state, in which the salinity in nodes on the fresh-salt water interface is constant with time, and the water flux in and out of the model is equal. The thing is, when I examine the model in terms of rate budget, I get huge volumes of water entering and leaving the system through the no-flow boundaries at the top and bottom of the model, as terms of distributed source/sink. Additionally, the y-component of the darcy flux isn't zero on those nodes where the boundary is set to no-flow. Both "source/sink" and "in/out transfer rate" in the material properties in the data panel are set to zero. The stranger thing is that even on nodes completely within the domain, far from the boundary, there are non-zero calculated values of rate budget. These nodes are on the fresh-salt interface, where density gradients are most significant. It is obviously related somehow to buoyancy term, but I can't figure out how is it possible that the system loses water from nodes away from the boundary or from nodes that are on a no-flow boundary.

Any help would be greatly appreciated.