For simplicity, we will write $f^2(x) = f(f(x))$, $f^3(x) = f(f(f(x)))$, etc.

First we show that $f$ must be injective. Suppose $f(z_1) = f(z_2)$. Replacing $z$ by $z_1$ or by $z_2$ in (1) leads to the same LHS expressions, hence the two RHS expressions must be equal as well, i.e. $z_1 + f(y) + f^2(x) = z_2 + f(y) + f^2(x)$ which implies $z_1 = z_2$.

We now make the RHS of (1) equal to $z + f^2(y) + f^2(x)$ in two different ways in order to benefit from the injectivity of $f$. Replacing $y$ by $f(y)$ gives $f(x + f^2(y) + f^2(z)) = z + f^2(y) + f^2(x)$. Replacing $x$ by $y$ and $y$ by $f(x)$ gives $f(y + f^2(x) + f^2(z)) = z + f^2(x) + f^2(y)$. Because the RHS is the same in both expressions, and $f$ is injective, it follows that
$$x + f^2(y) + f^2(z) = y + f^2(x) + f^2(z),$$
i.e. $f^2(x) - x = f^2(z) - z$ for all positive $x, z$.

This means that $f^2(x) - x = c$ is a constant that does not depend on $x$. Setting $t = 1 + f(1) + f^2(1)$, equation (1) with $x = y = z = 1$ implies $f(t) = t$, hence $f^2(t) = t$ and $f^2(t) - t = 0$. Therefore $c = 0$ and we have shown

$$f^2(x) = x \quad \text{for all } x > 0. \qquad (23)$$
Pick any $d > 0$ and put $x = d$, $y = f(d)$ and $z = d$ in (1), then use (23) to obtain $f(3d) = 3d$. This means that
$$f(x) = x \quad \text{for all } x > 0. \qquad (24)$$
This function clearly satisfies (1).