Solution (source code)

= Solution

Choose $f=x^2(1+y)+y^2z-1$. Since $f_z=y^2>0$ on the stated open domain, its upward <vector area element> is
$$
\boxed{d\mathbf S=\left(\frac{2x(1+y)}{y^2},\ \frac{x^2+2yz}{y^2},\ 1\right)dx\,dy,\qquad z=\frac{1-x^2(1+y)}{y^2}.}
$$
The graph becomes unbounded near parts of $y=0$; the flux integrals below are understood as limits from $y>0$.