Prove $\sqrt[2]{x+y}+\sqrt[3]{y+z}+\sqrt[4]{z+x} <4$

The condition $\;x^2+y^2+z^2+xyz=4\;$ can be rewritten as: $$ z^2+xy\cdot z + (x^2+y^2-4) = 0 \quad \Longrightarrow \quad z = - \frac{xy}{2} \pm \sqrt{\left(\frac{xy}{2}\right)^2-(x^2+y^2-4)} $$ From $z \ge 0$ it follows that: $$ z = - \frac{xy}{2} + \sqrt{\left(\frac{xy}{2}\right)^2-(x^2+y^2-4)} \quad \mbox{and} \quad x^2+y^2\le 4 $$ Therefore, after substitution of $\,z$ , we only have to investigate function behavior
$f(x,y) = \sqrt[2]{x+y}+\sqrt[3]{y+z}+\sqrt[4]{z+x} < 4$
in the area enclosed by x-axis, y-axis and a quarter of a circle, as depicted below.

Another proof without words is attempted by plotting a contour map of the function, as depicted. Levels (`nivo`) of these isolines are defined (in Delphi Pascal) as: 

nivo := min + g/grens*(max-min); { grens = 40 ; g = 1..grens }
The blackness of the isolines is proportional to the (positive) function values; they are almost white near the minimum and almost black near the maximum values. Maximum and minimum values of the function are observed to be: 
 2.54387743763872E+0000 < f < 3.91477606446737E+0000
The $\color{blue}{\mbox{blue}}$ spot is where $\,\left| f(x,y) - 4 \right| < 0.9$ . It is suggested by the rather large error $\,0.9\,$ that $\,4\,$ is not really the maximum. Indeed, upon refinement of the grid we find for the maximum numerically (double precision) a somewhat lower value: 
3.91477205860402 < 4