 Solution

Take the unit normal to point outwards. The divergence of the field isThe divergence theorem turns the closed surface integral into a volume integral over $[0,1{]}^3$:For a direct check of the surface integral, on $x=1$ the outward component is $F_x=2y$, contributing ${\int }_0^1{\int }_0^{1}2ydydz=1$. On $x=0$ it vanishes. On both $y=0$ and $y=1$, $F_y=xy(y-1{)}^2=0$. On $z=1$, $F_z=0$, whereas on $z=0$ the field component is $F_z=-1$ and the outward normal is $-\hat{\mathbf{z}}$, giving a contribution $1$. The six face contributions are therefore $1,0,0,0,0,1$, again totaling 2. The positive bottom-face contribution is a consequence of its outward orientation.