I think it is a rather hard problem to address. Non-uniformly sampled data are often interpolated into some equations, and a derivative of that is taken, as far as I know. There are many studies about numerical differentiation techniques in various situations. I found a past question that looks similar to your follow-up question, e.g. scicomp.stackexchange.com/q/480/23145

@Natasha It depends on how the sign of the boundary flux value is set. If the exiting boundary flux is given as negative, you need to flip the sign of the boundary flux equation as you suggested. In fact, it could make a better physical sense when it is looked as mass gain and loss of the domain. The sign of the second boundary condition equation in your original question, however, is not flipped.

It is amazing and unfortunate to see how inaccurate the information on the internet can be. I would rather recommend this Wikipedia document, en.wikipedia.org/wiki/…, to you. No matter where you look at, the equation should not vary.

The total fiux consists of two fluxes: The first one is advective flux that carries mass by the flow. This can be expressed in $vC$, which means the mass $C$ is transported along the direction of the flow. The second one is dispersive (or diffusive) flux that moves mass from higher concentration area to lower one. The direction of this flux is against the concentration gradient, so $-D \partial C / \partial x$.