When I was studying these topics, we had an approach based on the parametric representation of the surface rather than on differential forms. I will put an answer here since there still no other answer, but maybe it will not catch the whole generality.
Suppose that your surface is given by some parametrization $S = \{\mathbf r(u,v): u,v\in \Omega\}$ where $\Omega\subset \mathbb R^2$ and $(u,v)$ are parameters. Here $\mathbf r(u,v) = (x(u,v),y(u,v),z(u,v))$.
Then the surface integral of the first kind is given by $ \int\limits_S g(x,y,z)\,d S = \int\limits_\Omega g(\mathbf r(u,v))|\mathbf r_u\times\mathbf r_v|\,du \,dv $ where $\times$ denote the vector product in 3-dim vector space.
The unit normal to the surface $\mathbf n$ is given then as $ \mathbf n(u,v) = \frac{1}{|\mathbf r_u\times\mathbf r_v|}\left(\frac{\partial(y,z)}{\partial(u,v)},\frac{\partial(z,x)}{\partial(u,v)},\frac{\partial(x,y)}{\partial(u,v)}\right) $ where e.g. $\frac{\partial(y,z)}{\partial(u,v)} = \frac{\partial y}{\partial u}\frac{\partial z}{\partial v} - \frac{\partial y}{\partial v}\frac{\partial z}{\partial u}$ is a Jacobian. By the definition of the surface integral of the second kind as the vector flow through the surface, you get $ \int\limits_S f\,dy\,dz+g\,dz\,dx+h\,dx\,dy = \int\limits_\Omega (f,g,h)\cdot\left(\frac{\partial(y,z)}{\partial(u,v)},\frac{\partial(z,x)}{\partial(u,v)},\frac{\partial(x,y)}{\partial(u,v)}\right)\,du\,dv $ and you get the formula you are interested since $ \int\limits_S\mathbf F\cdot \mathbf n\,dS = \int\limits_\Omega\mathbf F\cdot\mathbf n|\mathbf r_u\times\mathbf r_v|du\,dv $