Evaluate $\iint_S x^2 dS$ where S is the paraboloid $z=x^2 + y^2$ for $0\le z \le 4$.
So far, I have parametrized my function as: $r(x,y) = \langle x, y, x^2 + y^2\rangle$ which gives me the partial sums: $$ r_x = \langle 1, 0, 2x\rangle$$ $$ r_y = \langle 0, 1, 2y\rangle$$ which gives me $ r_x \times r_y = \langle-2x, -2y, 1\rangle $. I'm a bit confused which orientation should I use and what my surface integral will be.
Rather than keeping the variables $x$ and $y$, consider parameterizing $S$ in cylindrical coordinates as
$$\vec s(r,\theta)=\langle r\cos\theta,r\sin\theta,r^2\rangle$$
which follows from the change of coordinates,
$$\begin{cases}x=r\cos\theta\\y=r\sin\theta\\z=x^2+y^2\end{cases}$$
where $0\le r\le2$ and $0\le\theta\le2\pi$. (Keeping the parameterization in terms of $x$ and $y$ would make the limits of integration more annoying to deal with.)
Now let the normal vector (again, direction doesn't matter) be
$$\vec n=\vec s_\theta \times \vec s_r=\langle 2r^2\cos t,2r^2 \sin t,-r\rangle$$
which makes the surface element
$$\mathrm dS=\|\vec n\|\,\mathrm dr\,\mathrm d\theta=r\sqrt{4r^2+1}\,\mathrm dr\,\mathrm d\theta$$
So the surface integral is
$$\begin{align} \iint_S x^2\,\mathrm dS&=\int_0^{2\pi}\int_0^2 (r\cos\theta)^2 r\sqrt{4r^2+1}\,\mathrm dr\,\mathrm d\theta\\[1ex] &=\left(\int_0^{2\pi}\cos^2\theta\,\mathrm d\theta\right)\left(\int_0^2 \, r^3 \sqrt{4r^2+1}\,\mathrm dr\right) \end{align}$$
which is easy to compute.
If you want to hold on to your original parameterization, you would have
$$\vec r(x,y)=\langle x,y,x^2+y^2\rangle$$
with $-2\le x\le2$ and $-\sqrt{4-x^2}\le y\le\sqrt{4-x^2}$. Then the normal vector would be as you found,
$$\vec n=\vec r_x\times\vec r_y=\langle-2x,-2y,1\rangle$$
so that the surface element is
$$\mathrm dS=\|\vec n\|\,\mathrm dx\,\mathrm dy=\sqrt{4x^2+4y^2+1}\,\mathrm dx\,\mathrm dy$$
and the integral would be
$$\iint_S x^2\,\mathrm dS=\int_{-2}^2\int_{-\sqrt{4-x^2}}^{\sqrt{4-x^2}}x^2\sqrt{4x^2+4y^2+1}\,\mathrm dy\,\mathrm dx$$