While tinkering with Mathematica I found that it could evaluate sums of the form $$\sum\limits_{k = 1}^{\infty} \frac{H_k H_{k+m}}{k^2}$$ for some nonnegative integer $m$. As it turns out (or the evidence is overwhelmingly in favor of it), this sum is always a linear combination of a constant term and the first three $\zeta$-values, i.e. $\zeta(2)$ through $\zeta(4)$. The coefficients of the $\zeta$-values were not too hard to figure out with the help of the OEIS and a bit of intuition. We have $$\sum\limits_{k = 1}^{\infty} \frac{H_k H_{k+m}}{k^2} = -C -\left(\sum\limits_{k = 1}^{m}\frac{1}{k^2}\right)\zeta(2) + 2\left(\sum\limits_{k = 1}^{m}\frac{1}{k}\right)\zeta(3) + \frac{17}{4}\zeta(4)$$ with just the non-zeta term $C$ being undetermined. For $m = 1,\ldots,9$ this constant takes on the values $$C = 0,\frac{1}{4},\frac{4}{9},\frac{341}{576},\frac{679}{960},\frac{25921}{32400},\frac{19879}{22680},\frac{95594629}{101606400},\frac{182134073}{182891520}$$ which... isn't all that insightful, to say the least. Visually, this looks like
Thanks to the comment by ho boon suan (and verification by Dr. Wolfgang Hintze) I was able to verify that indeed $$C = \frac{1}{2} \sum\limits_{j = 1}^{m} \frac{H_{j-1}^2+H_{j-1}^{(2)}}{j^2}$$
and thus
$$\sum\limits_{k = 1}^{\infty} \frac{H_k H_{k+m}}{k^2} = -\sum\limits_{j = 1}^{m} \underbrace{\frac{H_{j-1}^2+H_{j-1}^{(2)}}{2j^2}}_{=\frac{1}{j^2(j-1)!} \begin{bmatrix}j \\ 3\end{bmatrix}} -H_m^{(2)}\zeta(2) + 2H_m\zeta(3) + \frac{17}{4}\zeta(4)$$
As a side product of an auxiliary calculation, we also have the identity
$$\sum\limits_{k = 1}^{\infty} \frac{H_k}{k^2(k+j)} = \frac{H_{j-1}^2+H_{j-1}^{(2)}}{j^2} - \frac{1}{j^2} \zeta(2) + \frac{2}{j}\zeta(3)$$
for arbitrary $j \in \mathbb{N}$.