\begin{align*}
S &= \sum_{k \geq 1} \frac{1}{k^2} \sum_{r>s\geq 1} \frac{1}{(r^2 + s^2)^2} \\
&= \frac{\pi^2}{6} \sum_{r>s\geq 1} \frac{1}{(r^2 + s^2)^2}
\end{align*}
Note that $\frac{1}{(r^2 + s^2)^2}$ is unchanged by $\{r \mapsto -r\}$, $\{s \mapsto -s\}$, and $\{r \leftrightarrow s\}$. Applying these symmetries, we obtain eight copies of the region $r > s \geq 1$ in the $r$-$s$ plane. It lacks the four diagonals,
each with value $\displaystyle \sum_{s \geq 1} \frac{1}{(s^2 + s^2)^2} = \frac{\pi^4}{360}$, the four coordinate rays starting at $\pm 1$, each with value $\displaystyle \sum_{s \geq 1} \frac{1}{(0^2 + s^2)^2} = \frac{\pi^4}{90}$, and the origin, which we will continue to exclude. Placing a prime on the summation sign to indicate that we omit the origin,
$$ \sideset{}{'}\sum_{r,s = -\infty}^\infty \frac{1}{(r^2 + s^2)^2} = 8 \sum_{r > s \geq 1} \frac{1}{(r^2 + s^2)^2} + 4 \cdot \frac{\pi^4}{360} + 4 \cdot \frac{\pi^4}{90} $$
Borwein and Borwein [1] (p. 291) (see also (38) at Double Series on MathWorld) show
$$ \sideset{}{'}\sum_{r,s = -\infty}^\infty \frac{1}{(r^2 + s^2)^2} = 4 \beta(2) \zeta(2) = \frac{2}{3} \pi^2 K \text{,} $$
where $\beta$ is Dirichlet's Beta function, $\zeta$ is the Riemann Zeta function, and $K$ is Catalan's constant.
Utilizing these facts,
$$ S = \frac{\pi^2}{6} \cdot \frac{12 K \pi^2 - \pi^4}{144} = 0.1264945807\dots \text{.} $$
[1] Borwein, Jonathan M. and Peter B. Borwein, "Pi and the AGM: A Study in Analytic Number Theory and Computational Complexity", New York: Wiley, 1987.
P.S. Mathematica 11.3 does not evaluate the triple sum symbolically. It numerically estimates $0.126494$, which is comfortably close to the above result.
