According to [this answer][1],
$$z_\infty:=\lim_n z_n=I:=\int_0^\infty F(s)G(s)\,ds,$$
where 
$$F(s):=\prod_{k=1}^\infty\frac{1}{\sqrt{1+2s/k^3}},\quad G(s):=\sum_{k=1}^\infty\frac k{k^3+2s}.$$

Mathematica can express $F$ and $G$ in terms of functions built-in in Mathematica, and then the Mathematica command NIntegrate numerically evaluates $z_\infty=I$ as $\approx1.99218$ -- close to $2$, but not $2$; see the image of the corresponding Mathematica notebook below. 

Using the facts that (i) $F$ and $G$ are positive, decreasing, and convex, and hence $FG$ is so, and that (ii) Mathematica can find the values of all its built-in functions with any degree of accuracy, it should be rather straightforward to show that $I<2$. 

---

[![enter image description here][2]][2]

---

  [1]: https://math.stackexchange.com/a/4228443/96609
  [2]: https://i.sstatic.net/cIR4N.png