Let $\sigma_0(n)$ be the divisor counting function $$\sigma_0(n) = \sum_{d \vert n} 1.$$ I'm interested in the convolution sum $$ S(n) := \sum_{k=1}^{n-1} \sigma_0(k) \sigma_0(n-k)$$ I ran some quick numerical experiments and it looks like for odd primes $p$, the convolution sum $S(n)$ is equidistributed* mod $p$: $$\lim_{X \to \infty} \dfrac{\lvert \{ n<X: S(n) \equiv a \mod p \} \rvert}{X} = \dfrac{1}{p}.$$ Is this true? And if so, could anyone sketch a proof / provide a reference?
*This is a bit of a lie. In fact, the data seems to suggest that certain residue classes occur more often than others, see this question. But I've stated my observation in this (incorrect) way to avoid making the question too long. In any case, I'd like to learn of any techniques I could use to get a handle of $S(n)$ mod $p$.