The surface, z = x^2 + y^2, can be written in terms of parameters,u, v, as
x= u, y= v, z= u^2+ v^2 so that any point on the surface can be written in terms of the position vector \vec{r(u, v)}= u\vec{ix}+ v\vec{j}+ (u^2+ u^2)\vec{k}.
The derivatives with respect to the parameters are
\vec{r_u}= \vec{i}+ 2u\vec{j} and \vec{r_v}= \vec{j}+ 2v\vec{k}
Those are two vectors tangent to the surface, and their cross product,
\vec{r_v}\times\vec{r_u}= 2u\vec{i}+ 2v\vec{j}+ \vec{k}
times dudv, is the "vector differential of surface area". Its length, dS= \sqrt{4u^2+ 4v^2+ 1}dudv is the "differential of surface area".
So \int |xyz|dS= \int\int |uv(u^2+ v^2)(\sqrt{4u^2+ 4v^2+ 1}dudv
Seeing all those sum of squares I would be inclined to put it in "polar coordinates with u= rcos(\theta), v= r sin(\theta). z goes from 0 up to z= u^2+ v^2= r^2= 2 so that becomes
\int_{r= 0}^\sqrt{2}\int_{\theta= 0}^{2\pi} |r^4cos(\theta)sin(\theta)|\sqrt{4r^2+ 1}r drd\theta
To handle the absolute value, you will need to determine where cos(\theta)sin(\theta) is positive and where negative.
